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This  final  report  on  the  formulation  of  EPIC -3,  a 
three-dimensional  computer  code  for  dynamic  analyses 
of  high-velocity  impact  problems,  was  prepared  by 
Honeywell  Inc. , Defense  Systems  Division,  for  the 
U.S.  Army  Ballistics  Research  Laboratories,  Terminal 
Ballistics  Laboratory,  under  contract  DAAD05-76-C-0747, 
The  period  covered  by  the  report  is  January  1976  to 
February  1977, 

The  author,  Gordon  R.  Johnson,  would  like  to  thank 
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SECTION  I 

INTRODUCTION  AND  SUMMARY 


This  report  documents  a three-dimensional  computer  code  which  can  be 
used  for  dynamic  analyses  of  high-velocity  impact  problems.  The  code, 
EPIC- 3 (Elastic-Plastic  Impact  Calculations  in  ^ dimensions),  is  based 
on  a Lagrangian  finite -element  formulation  where  the  equations  of  motion 
are  integrated  directly  rather  than  through  the  traditional  stiffness  matrix 
approach.  Nonlinear  material  strength  and  compressibility  effects  are 
included  to  account  for  elastic-plastic  flow  and  wave  propagation.  Although 
this  code  is  arranged  to  provide  solutions  for  projectile-target  impact 
problems,  the  basic  formulation  is  valid  for  a wide  range  of  problems 
involving  dynamic  responses  of  continuous  media. 

The  EPIC -3  code  has  material  descriptions  which  include  strain  hardening, 
strain  rate  effects,  thermal  softening  and  fracture.  Geometry  generators 
are  included  to  generate  quickly  flat  plates  and  rods  with  blunt,  rounded  or 
conical  nose  shapes.  It  has  the  capacity  to  include  multiple  sliding  surfaces, 
it  is  restartable,  and  it  provides  cross-sectional  plots  of  the  deformed 
geometry, 

A desirable  characteristic  of  this  technique  is  that  there  is  no  need  to  pro- 
vide an  orderly  arrangement  of  nodes  as  is  required  for  finite -difference 
techniques.  Complex  geometrical  shapes  can  be  represented  simply  by 
providing  an  adequate  assemblage  of  elements  to  represent  the  geometry. 
Various  boundary  conditions  can  also  be  represented  in  a straightforward 
manner.  Another  desirable  characteristic  is  that  this  technique  is  for- 
mulated for  a tetrahedron  element,  which  is  well  suited  to  represent  the 
severe  distortions  which  often  occur  during  high-velocity  impact.  A 
tetrahedron  element  also  provides  a state  of  constant  strain  such  that  all 


1 


material  in  an  element  behaves  uniformly.  This  allows  for  an  accurate  and 
convenient  selection  of  a constant  stress  within  the  element. 

This  report  includes  a complete  description  and  formulation  of  the  EPIC -3 
computer  code.  Detailed  instructions  for  using  this  code  are  also  provided. 
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SECTION  II 
FORMULATION 


The  EPIC“3  computational  technique  is  shown  schematically  in  Figure  1. 
The  first  step  in  the  process  is  to  represent  the  geometry  with  tetrahedron 
elements  having  specific  material  characteristics.  Then  the  distributed 
mass  is  lumped  at  the  nodes  (element  corners),  and  initial  velocities  are 
assigned  to  represent  the  motion  at  impact. 
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Bhgure  1.  Basic  Computational  Technique 


After  the  initial  conditions  are  established,  the  integration  loop  begins  as 
shown  in  Figure  1.  The  first  step  is  to  obtain  displacements  and  velocities 
of  the  nodes.  If  it  is  assumed  the  lines  connecting  the  nodes  (element  edges) 
remain  straight,  then  the  displacements  and  velocities  within  the  elements 
must  vary  linearly.  From  these  displacements  and  velocities,  the  strains 
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and  strain  rates  within  the  elements  can  be  obtained.  Since  the  strains  and 
strain  rates  are  derivatives  of  linear  displacement  and  velocity  functions 
within  the  elements,  the  resulting  strains  and  strain  rates  are  constant 
within  the  elements. 

The  stresses  in  the  elements  are  determined  from,  the  strains,  strain  rates, 
internal  energies  and  material  properties.  Since  the  strains  and  strain  rates 
are  constant  within  the  elements,  the  stresses  are  also  constant.  If  the 
strains  are  in  the  elastic  range,  the  elastic  stresses  can  be  obtained  in  a, 
straightforward  manner  from  Hooke’s  law.  For  larger  strains,  involving 
plastic  flow,  the  stresses  are  obtained  by  combining  plastic  deviator  stresses 
with  hydrostatic  pressure.  The  plastic  deviator  stresses  represent  the  shear- 
strength  capability  of  the  material,  and  the  hydrostatic  pressure  is  obtained 
from  the  volumetric  strain  and  internal  energy  of  the  element.  An  artificial 
viscosity  is  also  included  to  damp  out  localized  oscillations  caused  by  repre- 
senting continuous  media  with  lumped  masses. 

After  the  element  stresses  are  determined,  it  is  necessary  to  obtain  concen- 
trated forces  at  the  nodes.  These  forces  are  statically  equivalent  to  the  dis- 
tributed stresses  within  the  element  and  are  dependent  on  the  element  geometry 
and  the  magnitude  of  the  stresses.  When  the  concentrated  forces  are  applied 
to  the  concentrated  masses,  the  nodal  accelerations  are  defined,  and  the 
equations  of  motion  are  applied  to  determine  new  displacements  and  velocities. 
The  integration  loop  is  then  repeated  until  the  time  of  interest  has  elapsed. 

An  important  additional  feature  of  the  basic  technique  is  the  ability  to  repre- 
sent sliding  between  two  surfaces.  This  is  accomplished  with  a momentum 
exchange  principle  which  allows  for  closing  and  separation  of  the  two  surfaces,. 
It  should  be  noted  that  the  integration  time  increment  must  be  properly  con- 
trolled to  prevent  numerical  instability.  This  is  accomplished  by  limiting 
the  time  increment  to  a fraction  of  the  time  required  to  travel  across  the 
minimum  altitu.de  of  the  element  at  the  sound  velocity  of  the  material.  This 
also  ensures  that  the  time  increment  is  less  than  the  lowest  period  of  vibra- 
tion of  the  system. 
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2.1  GEOMETRY 


A typical  tetrahedron  element  is  shown  in  PTgure  2.  It  is  geometrically 
defined  by  nodes  i,  j,  m,  and  p.  The  formulation  is  based  on  nodes  i,  j, 
and  m being  positioned  in  a counterclockwise  manner  when  viewed  from 
node  p.  For  the  initial  geometry,  the  coordinates  of  node  i are  defined  by 
x^,  V?,  and  z?.  The  coordinates  of  node  i for  the  displaced  geometry  are 
given  by  x.,  y.,  and  z..  In  all  cases,  the  superscript  "o"  is  used  to  designate 
the  initial  condition.  The  displacement  of  node  i along  the  x axis  is  u^  = x^ 

O 1 C J ^ 

-x^,  and  the  y and  z displacements  are  given  by  v-  y.  - y-  and  w.  - z,^  - z^  . 
The  lumped  mass  at  each  of  the  four  nodes  is 

M.'  = = M'  =M'=-i-V°p°  (1) 

I ] m p 4 ^ 

where  and  are  the  initial  volume  and  density  of  the  element.  The 
individual  elements  can  be  assembled  by  use  of  geometry  generators  des- 
cribed in  subsection  4.  1,  When  all  elements  are  completely  assembled  to 
form  a continuous  medium,  the  total  mass  at  node  i is 


2.2  STRAINS  AND  STRAIN  RATES 

The  strains  are  obtained  from  the  initial  geometry  of  the  element  and  the 
displacements  of  the  nodes.  If  it  is  assumed  the  displacements  vary  linearly 
between  the  nodes,  the  x,  y and  z displacements  (u,  v,  w)  can  be  expressed 
as 


u = 

“l 

(2) 

V = 

“5 

+ QTgX  + arjY  + cZgZ 

(3) 

w = 

“^9 

+ +^117  +0-12Z 

(4) 
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Figure  2,  Geometric  Properties  of  a Tetrahedron  Element 


where  a.  . , - are  geometry-  and  displacement- dependent  constants. 

It  is  possible  to  solve  for  ' ' ' ^4  substituting  the  x displacements 
and  initial  coordinates  of  the  four  nodes  into  Equation  (2),  This  gives  four 
equations  and  four  unknowns  so  that  the  constants  , a^)  can  be  evaluated. 

Equation  (2)  can  then  be  expressed  in  terms  of  element  geometry  and  nodal 
displacements. 
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where  the  initial  volume 

is 
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and  the  other  geometry- dependent  constants  are 

o 

X. 

3 

o 

y- 

o 

z. 

3 

0 

a.  = 

1 

o 

X 

m 

o 

V 

m 

o 

z 

m 

(6b) 

o 

X 

p 

o 

V 

-p 

o 

z 

p 

1 

o 

y- 

3 

o 
z . 
3 

ii 

o 

XI 

1 

o 

y 

m 

o 

z 

m 

(6c) 

1 

0 

y 

o 

z 

P 

1 

o 

X. 

3 

o 
z . 

3 

II 

o 

o 

1 

o 

X 

m 

o 

z 

m 

(6d)  ’ 

1 

o 

X 

p 

0 

z 

P 

d°=. 


o 

o 

x. 

y-i 

3 

j 

o 

o 

x 

y 

m 

m 

o 

o 

X 


y. 


(6e) 


The  remaining  geometry- dependent  constants  (a^^,  b?,  d?,  etc,  ) are 

obtained  by  a systematic  interchange  of  signs  and  subscripts.  The  y and  z 
displacements  in  Equations  (3)  and  (4)  are  obtained  in  a similar  manner  and 
are  identical  to  Equation  (5)  except  the  x displacements  are  replaced  by  the 
y and  z displacements. 


After  the  displacements  are  obtained,  it  is  possible  to  determine  the  strains 
within  the  elements  from  the  following  equations  (Ref.  1): 


€ 

X 


e 
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e 

z 


9u  , 1 

3x  ■ 2 


ivi  _L 
ay  2 
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'Bu 

w ^ . 

/- 

dw 

iBxj 
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Bu  ^ 

^y 

' 

Bw 

^y) 

^y) 

^ + 1 
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az| 


S z 
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(V) 

(8) 

(9) 


5u  ^ ^ 5u5u  ^ dw^W 

Sy  'dx  SxSy  SxSy  Sx3y 


(10) 


xz 


5u  ^ 5w  ^u^u  _j_  'dy'dy  ^ 5w^w 

c)z  3x  3xSz  BxSz  3xSz 


(11) 


3y  dw  5udu  dySy  dw5w 
Bz  By  ByBz  ByBz  ByBz 


(12) 


e 
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(13) 
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The  normal  strains  are  represented  by  Equations  (7)  through  (9),  the  shear 
strains  by  Equations  (10)  through  (12),  and  the  volumetric  strain  by  Equation 
(13).  The  higher-order  terms  are  required  to  represent  the  effect  of  element 
rotation. 


The  strains  in  Equations  (7)  through  (12)  are  determined  as  functions  of  the 
displacements  by  taking  partial  derivatives  of  Equations  (2)  through  (4), 
and  substituting  them  into  the  strain  equations.  These  strains  consist  of 
derivatives  of  linear  equations  and  are  therefore  constant  (functions  of  . 

o^ly)  within  the  element.  It  is  also  necessary  to  use  an  equivalent  strain 
which  is  expressed  as 


e = 


(e 


X 


e )^  + (e 

y X 


e )^  + (e  - e +|-(y  ^ + y ' 

z y z 2 xy  xz 


+ y ) 

yz 


(14) 


This  strain  provides  an  overall  measure  of  the  distortion  of  the  element  and 
can  be  used  to  obtain  strain  hardening  stresses  which  occur  during  plastic 
flow  (Ref,  2). 


The  strain  rates  are  obtained  in  a manner  similar  to  that  used  to  determine 
the  strains.  Since  the  velocities  (u,  v,  w)  also  vary  linearly  between  the  nodes. 


the  basic  form  of  Equations_  (2 ) through  (12)  is  used,  except  the  nodal  veloci- 
ties are  substituted  for  the  nodal  displacements  and  the  higher-order  terms 
are  dropped.  Also,  the  displaced  geometry  is  used  so  that  the  geometry- 
dependent  constants  in  Equation  (6)  are  obtained  from  the  displaced  coordinates 
(Xp  y*.  z.,  etc.)  rather  than  the  initial  coordinates  (x?,  y^,  z°,  etc,).  The 
resulting  constants  for  the  displaced  geometry  are  designated  as  a.,  b.,  c. , 
d.,  etc.  The  normal  strain  rates  (e 
strain  rates  (e  , e 


y , y 

xy  'xz 


and  y 


, e , e ) are  adjusted  into  deviator 
, X y z 

e ),  and  the  shear  strain  rates  are  represented  by 


y^ 
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2,3  STRESSES 

The  stresses  in  the  elements  are  determined  from  the  strains,  strain  rates, 
internal  energies,  and  material  properties.  The  elastic  stresses  can  be 
obtained  directly  from  the  strains  by  use  of  Hooke's  Law  (Ref,  2): 


(15) 

(16) 

(17) 

(18) 

(19) 

(20) 


'mal  stresses,  and  t , t and  t the  shear 

xy  xz  yz 

stresses.  Lame’s  elastic  constants  are  \ and  G,  and  the  artificial  viscosity, 
Q,  will  be  defined  later.  The  stresses  are  combined  to  form  an  equivalent 
stress  given  by 
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= ye  + 2GG  - 
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V X 
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y 

V y 
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= ye  -f  2Ge  - 
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V z 
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= Gy 
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xz 
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= Gy  „ 

yz 

yz 

where  d , o 

and  a are  the  i 

X y 

z 

~ a + (a  - a )^  + (a  - a )^  + 6 (T  ^ 
y X z ^ y z'  ' xy 


+ T 


xz 


+ T 


j 


(21) 


This  stress  is  analogous  to  the  equivalent  strain  of  Equation  (14)  inasmuch  as 
it  represents  an  overall  state  of  stress  within  the  element  (Ref.  2). If  a is  less  than 
the  tensile  yield  strength  of  the  material,  the  stresses  are  in  the  elastic 
range.  It  should  be  noted  that  the  elastic  stresses  are  dependent  on  strains 
which  are  relative  to  the  rotated  orientation  of  the  element.  Generally,  the 
element  rotation  is  very  small  while  it  is  in  the  elastic  range  and  no  com- 
pensation is  required.  If  significant  rotations  are  experienced  in  the  elastic 
range,  however,  a transformation  of  the  stresses  would  be  necessary. 
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Plastic  flow  begins  when  the  elastic  strength  of  the  material  is  exceeded.  In 
conjunction  with  the  plastic  flow,  allowance  is  also  made  to  include  the  effect 
of  viscosity.  For  these  cases  the  normal  stresses  are  obtained  by  combining 
the  plastic  and  viscous  deviator  stresses  with  the  hydrostatic  pressure  and 
artificial  viscosity.  This  gives 


- (P  + Q) 


a = s - (P  + Q) 

y y 


- (P  + Q) 


where  s , s and  s are  the  deviator  stresses,  P is  the  hydrostatic  pressure 
X y z 

and  Q is  the  artificial  vi,scosity. 


The  de viator  stresses  represent  the  shear  strength  characteristics  of  the 
material.  By  using  the  von  Mises  incremental  theory  of  plasticity  (Ref.  2) 
and  the  Navier-Strokes  equations  (Ref.  3).  the  deviator  stresses  (s^,  s^^,  s^) 
and  shear  stresses  (t  > given  by 
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The  first  terms  in  Equations  (25)  through  (30)  represent  the  plastic  strength 
where  S is  the  equivalent  tensile  strength  of  the  material  and  T is  the  equiva- 
lent strain  rate*  which  is  analogous  to  the  equivalent  strain  in  Equation  (14). 
The  second  terms  represent  the  viscous  effects  where  V is  the  absolute 
viscosity  oi;  the  inateriaL  For  the  case  of  no  viscosity  (V  = 0),  if  the  stresses 
in  Equations  (22)  through  (24)  and  Equations  (2  8)  through  (30)  are  substituted 
into  Equation  (21),  the  result  is  always  a " S.  It  should  be  noted  that  the  smm 
of  the  de viator  strain  rates  is  equal  to  zero,  and,  since  the  de viator  stresses 
are  directly  proportional  to  the  de  viator  strain  rates,  the  sum  of  the  de  viator 
stresses  is  also  equal  to  zero.  The  net  pressure  in  the  element  is  therefore 
not  affected  by  the  de  viator  stresses. 

The  equivalent  tensile  strength  of  the  material  may  be  dependent  on  many 
factors,  including  strain,  strain  rate,  pressure  and  temperature.  It  is  well- 
known  that  many  materials  behave  differently  under  dynamic  impact  than 
under  static  testing  conditions.  With  few  exceptions,  however,  a precise 
definition  of  material  behavior  under  dynamic  conditions  is  not  available. 

There  is  also  much  to  be  learned  about  fracture  characteristics  under  these 
dynamic  conditions.  This  version  of  EPIC-3  allows  the  equivalent  tensile 
strength  to  be  determined,  from 


S = S_  [l  + Cj  log  (e)]  [l  + C2H.  + CgH^]  [c^  + Cj-t; 


(31) 


In  this  formulation,  S"  is  generally  taken  to  be  the  static  stress,  which  is 
dependent  on  the  equivalent  strain  of  Equation  (14).  It  is  generally  necessary 
to  calculate  e on  a path* dependent  basis,  rather  than  the  path- independent 
basis  of  Equation  (14).  For  most  impact  problems,  however^  when  distortions 
are  severe  there  is  little  recovery  of  strains  and  this  is  an  adequate 
approximation. 
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The  three  bracketed  terms  allow  the  static  stress  be  altered,  based  on 
strain  rate,  pressure  and  temperature.  If  the  constants  are  = C2  ==  = 

= 0 and  C,  - 1,  0,  then  S = S- , which  is  the  strain  dependent  static  stress. 

5 4 ’ e ^ 

The  first  bracketed  term  can  increase  the  strengtii  due  to  the  equivalent 
strain  rate,  t,  described  earlier.  In  the  second  bracketed  term,  [jl  = p/p^  -1  = 
V‘^/V  -1,  Since  the  hydrostatic  pressure  will  be  shown  to  be  primarily  de- 
pendent on  p,  the  effect  of  this  term  is  to  increase  the  strength  due  to  pres- 
sure. The  final  bracketed  term  includes  the  temperat’.ire,  T,  of  the  element, 
which  is  obtained  from  the  nonrecoverabie  work  done  in  the  element  by  the 
artificial  viscosity  and  the  plastic  and  viscous  deviator  stresses.  This  term 
tends  to  reduce  the  strength.  The  formulations  for  the  nonrecoverabie  work 
and  temperature  are  presented  later  in  subsection  2.  7, 

Material  fracture  is  currently  dependent  on  the  equivalent  strain,  g , and  the 
volumetric  strain,  When  the  fracture  criterioan  has  been,  laiet,  the  equiva- 
lent tensile  stress  is  set  equal  to  zero,  so  that  no  shear  stresses  can  be 
developed  in  the  failed  element.  Likewise,  no  tensile  stresses  are  allowed 
to  develop.  The  net  result  is  that  a failed  element  tends  to  act  like  a liquid 
inasmuch,  as  it  can  develop  only  hydrostatic  compression  with  no  shear  or 
tensile  stresses.  Another  option  is  also  available  which  sets  all  stresses 
(including  the  pressure)  equal  to  zero. 

It  is  recognized  that  other  more  sophisticated  n.i.aterial  descriptions  Si.re  being 
developed.  It  should  be  noted  that  the  EPIC-3  formulation  is  ideally  suited 
to  incorporate  additional  material  models.  The  entire  tetrahedron  element 
is  in  a uniform  state  of  constant  strain  and  has  as Si> dated  with  it  the  para- 
meters such,  as  strains,  strain  rates,  internal  energy  and  temperatu.re , 

From  these  parameters,  the  stresses  in  the  ele.ment  can  readily  be  obtained 
by  various  selected  formulations. 

The  hydrostatic  pressure  is  dependent  on  the  volui.'.Letrii;  strain  and  the 
internal  ene rgy  in  the  element.  The  EPIC -3  code  uses  the  Mie-Gruneisen 
equation  of  state  in  the  following  form  (Ref,  4)' 
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(32) 


P = (KjH,  +K.2M.'"  +Kg|x^)  jl  - Jij  + (1  +^) 

The  internal  energy  per  original  unit  volume,  E , is  obtained  from  the 

S 

element  by  the  various  stresses  and  is  formulated  later  in  2.  7;  K, , EC  and 
Kg  are  material-dependent  constants,  and  r is  the  Grflneisen  coefficient. 

The  artificial  viscosity  is  combined  with  the  normal  stresses  to  damp  out 
localized  oscillations  of  the  concentrated  masses.  It  tends  to  eliminate 
spurious  oscillations  which  would  otherwise  occur  for  wave  propagation 
problems.  This  technique  was  originally  proposed  by  Von  Neumann  and 
Bdchtrnyer  (Ref.  5)  and  has  been  expanded  for  use  in  various  computer  codes 
(Ref,  6,  7).  It  is  expressed  in  terms  of  linear  and  quadratic  components 
and  is  applied  only  when  the  volumetric  strain  rate  is  negative: 

Q = C.  pc  h |e  I f Cmh^(0  for  e <0  (33) 

o V U Y V 

Q = 0 for  € ^ 0 

V 


where  is  the  sound  velocity  of  the  material  and  h is  the  minimum  altitude 
of  the  tetrahedron.  Typical  values  used  for  the  dimensionless  coefficients 

9 

are  C,  = 0.  5 and  ~ 4.0. 

L o 


The  minimum,  altitude  of  the  tetrahedron  can  readily  be  obtained  from  the 
previously  calculated  geometry- dependent  constants.  The  altitude  from  node 
i to  the  plane,  defined  by  the  other  three  nodes  of  the  tetrahedron  is 


6V 


2 

c.r  + d:' 

I I 


(34) 


The  other  altitudes  (li. , ly  , h ) are  obtained  by  appropriately  changing  the 

3 m p 

subscripts  of  the  geometry  dependent  constants. 
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The  sound  velocity  is  dependent  on  the  state  of  the  material  (elastic  or 
plastic)  and  is  defined  as  follows: 


c -\{X+2G)/p^^ 
s 


0.(K^  + + 3Kj|J.2)/p° 


(35a) 

(35b) 


2.4  CONCENTRATED  FORCES 


After  the  element  stresses  are  determined,  the  concentrated  nodal  forces 
can  be  obtained.  These  forces  are  statically  equivalent  to  the  distributed 
stresses  within  the  element  and  are  dependent  on  the  displaced  element 
geometry  and  the  magnitude  of  the  stresses.  The  forces  in  the  x,  y and  z 
directions  at  node  i of  an  element,  are  given  by 


= - 4 (b-CT  + c-T  + d.  T ) 
L 6 l X L XV  L XZ 


F?^  = -w(c.a  +b.T  +d.T  ) 
L 6 L y L xy  1 yz 


F^  = -y(d.a  +C.T  +b.T  ) 
L 6 1 z 1 yz  L XZ 


(36) 

(37) 

(38) 


The  geometry-dependent  constants  (b^,  d.)  are  again  identical  to  those 

used  for  calculation  of  the  strain  rates  and  altitudes.  The  forces  at  the  other 
nodes  are  readily  obtained  by  changing  subscripts.  The  net  forces  at  node  i 
2f?)  are  the  sum  of  the  node  i forces  from  each  element  which 
includes  that  node. 
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2.  5 EQUATIONS  OF  MOTION 

The  equations  of  motion  can  be  numerically  integrated  by  assuming  a con 
stant  velocity  for  each  time  increment.  The  acceleration  of  node  i in  the 
X direction  at  time  = t is 


..  t 

u- 


{39) 


The  new  constant  velocity  for  the  next  time  increment  is 


' t+  • t-  , « t . , 
u.  = u . + u . At 

i i 1 


(40) 


• t- 

where  u.  is  the  constant  velocity  for  the  previous  time  increment  and  At  is 
the  integration  time  increment.  Finally,  the  new  displacement  at  time  ^ t+At 
is 


t+At 

u. 

1 


i 


At 


(41) 


The  equations  of  motion  for  the  y and  z directions  have  a similar  form. 


It  should  be  noted  that  both  translational  and  rotational  momesita  are  conserved 
with  this  formulation.  The  concentrated  forces  of  Equations  (36)  through (3 8)  are 
in  static  equilibrium  so  that  the  net  forces  and  moments  associated  with  each 
element  are  equal  to  zero.  Therefore,  the  net  forces  and  moments  acting  on 
the  entire  system  are  also  equal  to  zero  if  there  are  no  external  restraints  or 
applied  forces.  Since  the  impulses  acting  on  the  nodes  are  simply  the  forces 
times  the  time  increment,  the  net  impulses  are  equal  to  zero,  and  the  net 
translational  and  rotational  momenta  of  the  system  are  not  altered. 
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The  integration  time  increment  used  for  the  equations  of  motion  is  given  by 


(42) 


2 2 

where  g = C Q/p,  h is  the  previously  determined  minimum  altitude,  and 
is  the  sound  velocity.  The  constant,  must  loe  less  than  unity  to  ensure 

that  At  is  always  less  than  the  time  required  to  travel  across  the  shortest 
dimension  of  the  tetrahedron  at  the  sound  velocity  of  the  material.  Use  of  a 
time  increment  significantly  larger  than  that  specified  by  Equation  (42)  will 
lead  to  numerical  instability  (Ref.  7).  The  EPIC -3  program  restricts  the 
time  increment  from  increasing  more  than  10  percent  per  cycle. 


2.  6 SLIDING  SURFACES 

It  is  sometimes  necessary  to  allow  for  sliding  to  occur  between  two  surfaces, 
A summary  of  the  Important  steps  for  the  sliding-surface  technique  follows: 

• Identify  a "master"  sliding  surface  defined  by  an  orderly 
arrangement  of  master  nodes, 

• Identify  "slave"  nodes  which  may  slide  along  the  master  surface, 

• For  each  integration  time  increment,  apply  the  equations  of 
motion  to  both  the  master  nodes  and  the  slave  nodes  in  the 
usual  manner, 

• For  each  slave  node,  find  the  triangular  plane  (defined  by 
three  master  nodes)  whose  pro;jection  contains  that  slave 
no  de , 

• Check  to  determine  if  there  is  interference  between  the 
slave  node  and  the  master  surface  (the  triangular  plane 
defined  by  three  master  nodes), 
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• If  there  is  interference,  place  the  slave  node  on  the  master 
surface  in  a direction  normal  to  the  appropriate  master 
triangular  plane, 

• Calculate  the  momentum  change  (impulse)  to  the  slave  node 
caused  by  placing  it  on  the  master  surface;  then  transfer 
this  lost  momentum  to  the  three  surrounding  master  nodes. 

The  first  step  in  the  process  is  to  define  the  master  surface.  It  is  geneirally 
desirable  to  define  this  surface  with  an  orderly  arrangement  of  nodes  since 
a systematic  search  is  required  on  this  surface.  An  example  of  this  is  shown 
in  the  upper  portion  of  Figure  3 where  a flat  plate  is  represented  by  a finite- 
element  model  of  nodes  and  elements.  The  z-axis  points  vertically  upward 
and  the  x and  y axes  are  in  a horizontal  plane,  fl.lthough  the  basic  principles 
involved  in  this  technique  are  independent  of  the  orientation  of  the  axes,  the 
search  technique  is  simplified  if  the  z coordinates  of  the  master  surface  are 
limited  to  a single-valued  function  of  x and  y (i,  e,  , any  vertical  line,  parallel 
to  the  z axis  must  not  pass  through  the  master  surface  at  more  than  one  point). 
Also,  for  this  discussion,  the  slave  surface  is  above  the  master  surface 
relative  to  the  z axis.  Due  to  the  symmetry  of  the  plate  about  a plane  con- 
taining the  X and  z axes,  only  half  the  plate  is  considered  as  shown.  The 
EPIC-3  code  uses  a master  surface  defined  by  M rows  of  nodes,  with  each, 
row  containing  N nodes.  Looking  downward  from  the  positive  z axis,  the 
second  row  of  nodes  from  N + 1 to  2N,  is  to  the  left  of  the  first  row  of  nodes 
going  from  1 to  N.  Likewise,  the  third  row  is  to  the  left  of  the  second  row, 
etc.  For  the  specific  case  shown  in  Figure  3,  N *-  2 3 and  M = 4,  for  a total 
of  92  master  nodes. 

It  is  also  necessary  to  identify  slave  nodes  for  the  other  sliding  surface. 

Since  these  nodes  do  not  require  an  orderly  arrangement,  it  is  usually 
convenient  to  designate  the  more  geometrically  com.plex  surface  as  the 
slave  surface.  It  is  also  desirable  to  I'estrict  the  slave  node  spacing  from, 
becoming  significantly  greater  than  the  master  node  spacing,  because  this 
could  introduce  localized  deformations  in  the  master  surface  at  the  slave 
node  locations.  Likewise,  the  mass  of  the  slave  nodes  should  not  be 
significantly  greater  tha.n  the  mass  of  the  .master  nodes,  since  this  could 
result  in  unrealistically  high  velocities  of  the  master  nodes  at  impact, 
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X-Y  PROJECTION  OF  DISPLACED  MASTER  SURFACE 

lire  3.  Ma.ster  Surface  Geometric  Representation 


For  each  integratiori  loop,  the  equations  of  JXi.otion  are  applied  to  both  the 

master  and  slave  nodes  in  the  usual  manner  described  in  subsection  2.  5. 

It  is  then  necessary  to  check  each  slave  node  to  determine  if  it  has  passed 
through  the  master  surface,  thus  causing  interference.  To  begin,  the 
extreme  values  of  the  entire  master  surface  are  determined;  if  the  slave 
node  is  not  within  this  region,  there  can  be  no  interference.  If,  for  instance , 
the  z coordinate  of  the  slave  node  is  greater  than  the  maximum  z coordinate 
of  the  master  surface,  there  can  be  no  interference,  and  the  check  for  that 
particular  slave  node  is  complete. 

If  the  slave  node  is  within  the  confines  of  the  master  surface,  a search 
through  the  master  surface  must  be  initiated.  It  is  for  this  search  that  it 
is  desirable  to  have  an  orderly  arrangement  of  master  nodes.  The  lower 
portion  of  Figure  3 shows  the  projection  of  an  arbitrarily  displaced  master 
surface  on  the  x-y  plane.  The  triangular  grid  is  consistent  with  the  triangular 
faces  of  the  tetrahedron  elements.  Since  z is  a single- valued  function, 
there  is  no  overlapping  of  triangles  on  the  x-y  projection.  The  objective  of 
the  search  is  to  determine  which  triangle  contains  the  x~y  projection  for 
each  slave  node.  The  search  begins  by  considering  the  triangles  between 
the  first  two  rows  of  nodes.  With  reference  to  the  lower  portion  of  Figure 
3,  the  first  triangle  is  defined  by  nodes  1,  2,  24,  the  second  by  nodes  2, 

25,  24,  and  the  third  by  nodes  2,  3,  25,  etc.  If  the  slave  node  projection 
is  not  found  by  the  last  triangle  in  the  row  (nodes  23,  46,  45),  the  pattern 
is  repeated  in  the  second  row  of  triangles  beginning  with  the  triangle  defined 
by  nodes  24,  25,  47, 


An  efficient  way  to  check  if  the  slave  node  projection  is  within  a specific 
triangle  is  to  determine  the  distance  from  the  slave  node  to  each  of  the  lines 
defining  the  three  sides  of  the  triangle.  The  distance  to  line  i-j,  in  the  x-y 
plane,  is  given  by 


3 


xy 


(43) 
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where  [3  - y.  - y.,  (3  = x.  - x.  and  /3  = x.y.  - x.y.-  The  coordinates  of 

1 1 2 3 1 3 r^3  j-i 

master  node  i are  x. , y.  and  z.,  and  the  coordinates  of  the  slave  node  are 
y\  and  z'  If  master  node  j is  counterclockwise  from  node  i,  when 

O O S 

looking  downward  from  the  positive  z axis,  and  5 ^ is  positive  for  all  three 
lines  of  the  triangle,  then  the  slave  node  is  contained  within  the  triangle. 
This  is  also  illustrated  in  Figure  3 where  the  slave  node  is  included  in  the 
triangle  defined  by  nodes  40,  41  and  63. 


After  the  appropriate  triangle  is  identified,  it  is  necessary  to  determine  if 
there  is  interference.  The  equation  of  a plane  through  nodes  i,  j and  k 
(taken  counterclockwise,  when  looking  downward  from  the  positive  z axis) 
is 

Ax  + By  + Cz  + D = 0 (44) 
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If  the  coordinates  of  the  slave  node  are  j,'  z'  the  corresponding  2: 

o S o 

coordinate  of  the  plane,  z , can  be  obtained  by  rearranging  Equation  (44): 

P 


z 

P 


D + Ax  ^ + By  ^ 
s s 

C 


(46) 


If  z (master)  is  less  than  (slave),  the  surfaces  are  separated  at  this 
p ' s 

slave  node  location,  and  the  final  coordinates  of  the  slave  node  (x  , , z^,) 

are  simply  equal  to  the  initial  coordinates  (x'  y^,  z ')  determined  from 

s s s 

the  equations  of  motion.  If,  however,  z is  greater  than  z'  there  is 

p b 

interference,  and  the  sliding  process  must  be  continued. 


If  there  is  interference,  the  next  step  is  to  place  the  slave  node  on  the 
master  surface  in  a direction  normal  to  the  master  surface.  This  tech- 
nique is  illustrated  in  Figure  4 where  a triangular  plane  is  defined  by 
master  nodes  i,  j and  k.  The  initial  slave  node  position  (x'  y'  z ')  lies 
below  the  triangular  plane,  as  shown  in  the  comparison  between  z^  and  z^. 
A vector  is  shown  which  passes  through  the  initial  slave  node  position  and 

is  normal  to  the  triangular  plane.  The  three  direction  cosines  of  this 

2 2 

vector  are  readily  obtained  by  dividing  A,  B and  C by  VA  + B + C , The 
slave  node  is  moved  a distance,  6^^,  along  this  normal  vector  until  it 
intersects  the  triangular  plane.  This  normal  distance  is  defined  by 


6 


N 


Ax'  + By'  + Cz'  + D 

S ''  s s 

~\Ja^  + 


(47) 


The  final  coordinates  of  the  slave  node  are  + Ax,  Yg  “ Jg  + 

z = z'  + Az,  where  Ax,  Ay  and  Az  are  obtained  by  multiplying  by  the 
s s 

appropriate  direction  cosines  of  the  normal  vector. 


After  geometric  compatibility  has  been  achieved  by  placing  the  slave  node 
on  the  master  surface,  it  is  necessary  to  adjust  the  equations  of  motion  to 
account  for  this  change.  Since  the  geometric  changes  are  made  only  in  the 
direction  normal  to  the  master  plane,  the  equations  of  motion  are  altered 
only  in  this  normal  direction.  No  changes  are  made  to  the  equations  of 
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Figure  4.  Slave  Node  Placement  onto  Master  Surface 
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motion  parallel  to  the  master  surface.  This  has  the  effect  of  assuming  a 
frictionless  surface. 

In  the  process  of  moving  the  slave  node  to  the  master  surface,  an  effective 
velocity  change  is  imposed  on  the  slave  node.  This  loss  of  velocity  in  the 
normal  direction  is  AV^  = 6j^/At,  where  At  is  the  integration  time  incremer 
for  the  previous  cycle.  This  velocity  change  also  gives  a momentum  change 
along  the  normal  direction  equal  to  AV^^,  where  is  the  mass  of  the 
slave  node.  This  momentum,  is  then  transferred  to  the  three  surrounding 
master  nodes.  If  and  Rj^  represent  the  fractions  of  the  momentum 

to  be  transferred  to  the  three  nodes,  a summation  of  the  translational  and 
rotational  momenta  of  the  three  nodes  gives  the  following  equations: 

Rj  + R,j  + = 1 (48) 

R.  (X.  - Xg)  + R.  (X.  - Xg)  + R^  (x^  - = 0 (49) 

R;  (>'i  - yg)  + H-  (y^.  - yg)  + \ (y^  - y^)  = o (50) 

These  three  equations  and  three  unknowns  define  the  appropriate  momentum 
distribution  to  the  three  nodes.  This  momentum  is  applied  to  the  master 
nodes  in  the  form  of  instantaneous  velocity  changes,  parallel  to  the  normal 
vector  but  in  the  opposite  direction.  For  instance,  the  normal  direction 
velocity  change  to  node  i is 

AV.  = R.  AV,,/M.  (51) 

1 1 s N 1 

where  M.  is  the  mass  of  master  node  i,  and  M is  the  mass  of  the  slave 

node.  As  was  done  for  the  displacements,  the  specific  velocity  changes 

in  the  x,  y and  z directions  are  obtained  by  multiplying  AV.  by  the  appro- 

1 

priate  direction  cosines. 

This  is  an  extremely  accurate  technique  because  translational  momenta  are 
absolutely  conserved  in  all  three  directions.  Small  errors  are  introduced 
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into  the  center  of  gravity  (eg)  positions  when  the  slave  node  is  moved  to  the 
master  surface,  since  there  is  no  corresponding  movement  of  the  master 
nodes,  which  receive  only  instantaneous  velocity  changes.  This  eg  error 
also  introduces  small  errors  into  the  rotational  momenta.  The  magnitude 
of  these  errors  is  identified  in  the  example  problem  of  subsection  2.  9, 

2,  7 SYSTEM  PARAMETERS 

It  is  often  desirable  to  examine  system  parameters  such  as  eg  positions, 
translational  and  rotational  momenta,  kinetic  energy,  plastic  work  and 
total  energy.  The  EPIC-3  code  calculates  these  parameters  for  the  pro- 
jectile, target  and  total  system  (projectile  and  target).  The  eg  position  in 
the  X direction  is  given  by 

X = 

where  M-  and  x^  are  the  mass  and  x coordinate  of  node  i.  The  summation 
includes  all  nodes  of  the  item  of  interest  (projectile,  target  or  total  system). 
The  eg  positions  for  the  y and  z axes  are  determined  in  a similar  manner, 

The  three  translational  momenta  (MV  , MV  , MV  ) are  obtained  in  a similar 

X y z 

manner,  with  the  momentum  in  the  x direction  given  by 

MV  -MM.u.  (53) 

where  uMs  the  x velocity  of  node  i and  the  summation  again  includes  all 
nodes  of  the  item  of  interest,  If  there  are  no  external  restraints  in  a 
particular  direction,  the  translational  momentum  in  that  direction  will 
remain  unchanged  for  the  total  system. 


EM.x. 
1 1 

S M . 

]. 


25 


The  three  rotational  momenta  (H  , H , H ) are  obtained  about  the  eg  of  the 

X y z 

item  of  interest.  The  rotational  momentum  about  the  x axis  is  defined  as 
H ^ I^M.[w.(y.  - y)  - v.  (z.  - z)]  (54) 

X X 1-  1 XX 

where  v-  and  z.  are  the  coordinates  of  node  i,  and  v.  and  w.  are  the  cor- 
1 11 

responding  velocities. 


The  kinetic  energy  is  next  determined  from 

KE  = y M . ( u^.  + V?  + ) (55) 

For  the  initial  conditions  it  is  assumed  there  is  no  internal  energy  which 
can  be  converted  to  kinetic  energy.  Therefore,  the  kinetic  energy  of  the 
system  cannot  increase  above  its  initial  level.  It  can,  however,  decrease 
by  being  converted  to  various  forms  of  internal  energy. 


The  eg  positions,  momenta  and  kinetic  energies  are  all  parameters  which 
are  determined  from  the  node  variables.  Various  forms  of  internal  energy 
are  associated  with  the  elements  and  are  required  for  the  determination  of 
the  temperature,  T,  in  Equation  (31)  and  the  hydrostatic  pressure,  P,  in 
Equation  (32).  The  internal  energy  in  an  element  is  approximately  deter- 
mined from 


E 


t +.^t 


= E*"+[se  +se  +se  +t  y 

^ -xy  -%.r  -XT  -XT  rr  rr  -vrTr'  ■ 


XX  y y z z xy'xy 


+ T y +T  y -(P+Q)0  ]At 
xz^xz  yz^jT. 


(56) 


where  E is  the  internal  energy  per  original  unit  volume  at  time  t+At. 

s 

The  deviator  stresses  (s  , s , s ),  shear  stresses  (t  , t , ),  hydro- 

X y z xy  xz  y z 

static  pressure,  and  artificial  viscosity  (P,Q)  are  those  obtained  at  time  t+At. 

Likewise,  the  strain  rates  (e  , e , e , y , y^^^,  e ) ^.re  for  this  same 

X y z xy  xz  yz  v 

time.  The  time  increment.  At,  is  for  the  previous  cycle  between  times  t and 
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t+At,  A more  precise  determination  of  the  internal  energy  could  be  made 
by  considering  the  stresses  and  strains  at  times  t and  t-At,  but  this  would 
require  a significant  increase  in  storage  requirenientvS, 

The  total  internal  energy  in  an  individual  element  at  time,  t,  is  where 

o ^ 

V is  the  initial  volum.e  of  the  element.  The  total  internal  energy  for  the 

system  is  the  sum  of  the  internal  energies  of  the  individual  elements.  The 

.total  system,  energy  is  the  sum  of  the  system  kinetic  energy  and  the  system 

internal  energy,  and  it  should  be  equal  to  the  initial  kinetic  energy.  A 

comparison  of  the  total  system  energy  ,with.  the  initial  kinetic  energy  gives 

an  indication  of  the  accuracy  of  the  approximate  formulation  of  Equation  (56), 

This  approximation  is  generally  very  accurate  for  high-velocity  impact 

problems  when  the  strains  are  constantly  increasing.  It  may  be  less  accurate, 

however,  for  problems  where  the  strains  are  not  constantly  increasing. 

The  effect  of  E^  in  Equation  (32)  can  be  eliminated  by  setting  the  Grtineisen 

coefficient  equal  to  zero. 


The  temperature  of  an  element  is  obtained  from  the  nonrecoverable  work 
done  on  the  element.  When  an  element  is  in  the  plastic  range,  the  work  per 
original  unit  volume,  E^,  is  identical  to  that  of  Equation  (56),  except  the 
effect  of  the  hydrostatic  pressure,  P,  is  not  included  since  its  v/ork  can  be 
recovered.  The  work  done  by  the  plastic  and  viscous  deviator  stresses, 
shear  stresses  and  artificial  viscosity  cannot  be  recovered  and,  therefore, 
causes  an  increase  in  temperature.  When  the  material  is  in  the  elastic 
range,  only  the  work  done  by  the  artificial  viscosity  cannot  be  recovered. 
The  temperature  of  an  element  is 


T - T 

o 


(5  7) 


where  T is  the  initial  temperature,  C is  the  specific  heat,  and  is  the 
o s 

initial  density.  The  temperature  only  occurs  in  Equation  (31),  and  its 
effect  can  be  eliminated  by  setting  - 1.0  and  C„  - 0. 


2.  8 SEVERE  DISTORTIONS 


It  was  previously  stated  that  a tetrahedron  element  is  well  suited  to 
represent  severe  distortions.  In  Reference  8 it  is  shown  that  a two- 
dimensional  triangular  element  formulation  is  much  better  suited  to 
represent  severe  distortions  than  is  a formulation  which  uses  a quad- 
rilateral grid.  This  is  due  to  a triangle's  resistance  to  allowing  a node  to 
cross  the  opposite  side  of  the  triangle  during  high  distortion.  It  is  apparent 
that  the  cross-sectional  area  of  a triangle  approaches  zero  as  a node 
approaches  the  opposite  side  of  that  triangle.  Since  the  hydrostatic  pressure 
is  inversely  proportional  to  the  volume  of  the  element,  as  the  volume 
approaches  zero,  the  pressure  approaches  infinity.  As  a result,  there  is 
a very  large  resistance  to  this  mode  of  distortion,  and  the  triangle  will  not 
go  through  a zero-volume  configuration,  provided  the  time  increment  is 
properly  controlled. 

The  quadrilateral  formulation  does  not  have  this  desirable  characteristic. 
Figure  5 shows  three  quadrilateral  configurations  which  have  equal  cross- 
sectional  areas,  A^,  Configuration  A shows  a rectangular  cross  section  with, 
nodes  i,  j,  rn  and  p defining  the  initial  geometry.  Configuration  B represents 
a deformed  element,  with  nodes  j and  m occupying  the  same  position.  A 
triangular  element  with  two  nodes  occupying  the  same  position  would  result 
in  a zero  cross-sectional  area.  For  the  quadrilateral  element,  however, 
nodes  i and  p can  be  displaced  to  allow  the  initial  area  to  be  retained.  Con- 
figuration C shows  overlapping  nodes,  and  the  net  area  of  the  quadrilateral 
is  equal  to  the  difference  of  the  two  triangular  cross  sections.  Even  this 
configuration  can  exist  with  the  initial  area.  When  this  type  of  distortion 
is  achieved,  however,  the  formulation  generally  breaks  down,  and  numerical 
instability  often  occurs.  The  deformations  represented  by  configurations 
B and  C usually  occur  when  the  grid  is  highly  distorted.  These  deforma- 
tions are  resisted  by  the  deviator  and  shear  stresses,  but  these  stresses 
are  limited  by  the  strength  of  the  material  and  do  not  provide  the  potential 
magnitude  of  resistance  as  does  the  hydrostatic  pressure  for  the  triangular 
elements. 
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The  three-dimensional  tetrahedron  element  used  in  the  EPIC- 3 code  is 
analagous  to  the  two-dimensional  triangular  element  with  respect  to  severe 
distortions.  The  lines  connecting  the  nodes  cannot  become  tangled  or 
overlapped  since  the  tetrahedron  must  also  go  through  a zero -volume  con- 
figuration (and  therefore  infinite  pressure)  before  this  can  occur.  This  is 
not  the  case  for  any  three-dimensional  element  with  more  than  four  sides. 
Since  many  impact  problems  produce  severe  distortions,  this  is  a very 
important  feature  of  the  tetrahedron  element  formulation. 


2.9  EXAMPLE 

An  example  problem  is  shown  in  Figure  6 to  illustrate  the  EPIC-3  com- 
putational technique.  It  involves  oblique  impact  (60  degrees  from  normal) 
of  a steel  sphere  onto  an  aluminum  plate  at  729  m/s.  The  sphere  diameter 
and  plate  thickness  are  both  equal  to  0.  635  cm.  The  yield  and  ultimate 
stresses  are  taken  to  be  the  static  values  of  2,  07  and  2,2  8 GPa  for  the 
steel  sphere,  and  0,  31  and  0.  45  GPa  for  the  aluminum  plate.  This  problem 
was  selected  because  test  data  in  Reference  9 indicate  the  sphere  should 
ricochet  off  the  surface  of  the  plate. 


The  initial  geometry  in  Figure  6 is  defined  by  831  nodes  and  3192  tetrahedron 
elements.  Note  that  the  geometry  of  the  sphere  can  be  adequately  represented 
with  the  tetrahedron  elements.  The  outer  nodes  of  the  sphere  are  the  slave 
nodes,  and  the  nodes  on  the  top  surface  of  the  plate  are  the  master  nodes  as 
shown  in  Figure  3,  The  deformed  shape  at  12  ps  after  impact  is  shown, 
and  this  is  the  time  at  which  the  vertical  component  of  the  sphere  velocity 
is  equal  to  zero.  Thereafter,  the  sphere  moves  upward,  away  from  the 
plate.  At  40  ps  after  impact,  the  sphere  has  separated  from  the  plate  and 
is  flying  freely.  This  problem  required  507  cycles  of  integration  and 
slightly  less  than  three  hours  of  central  processor  time  on  a Honeywell 
6080  computer. 
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The  cross-section  of  the  plate  at  40  [j-s  is  in  good  general  agreement  with 
the  test  data  reported  in  Reference  9.  It  is  anticipated  that  there  would  be 
even  better  agreement  if  the  finite-element  model  were  extended  in  the 
horizontal  directions  to  represent  a larger  plate,  if  a finer  network  of 
nodes  were  used,  and  if  dynamic  material  properties  were  used  in  place  of 
the  static  material  properties. 

It  was  previously  mentioned  that  slight  errors  in  rotational  momenta  and 
eg  positions  can  result  with  this  formulation.  Since  the  initial  path  of  the 
sphere  does  not  pass  through  the  eg  of  the  plate,  the  initial  system  (sphere 
and  plate)  has  a rotational  momentum  about  the  axis  normal  to  the  plane  of 
symmetry.  At  40  ps  after  impact,  when  the  sphere  is  flying  freely,  this 
rotational  momentum  is  within  1 percent  of  the  initial  value.  The  eg  errors 
associated  with  the  rigid-body  translation  of  the  system  are  significantly 
less  than  1 percent.  There  are  no  errors  in  the  translational  momenta 
since  these  parameters  are  conserved  in  the  basic  formulation. 

Additional  examples,  involving  wave  propagation,  severe  distortions  and 
material  fracture,  are  provided  in  References  10  and  11, 
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SECTION  III 

COMPUTER  PROGRAM  DESCRIPTION 


The  EPIC -3  computer  program  contains  the  formulation  presented  in 
Section  II.  A description  of  some  of  the  characteristics  of  the  computer 
program  is  given  in  the  following  subsections. 


3,  1 PROGRAM  ORGANIZATION  AND  SUBROUTINES 


The  organization  of  the  EPIC -3  code  is  shown  in  Figure  7,  It  consists  of 
the  main  program,  EPIC -3,  and  21  subroutines  for  a total  of  4344  cards. 

A plotting  program  is  also  available  to  give  two-dimensional  plots  of  the 
x-z  plane  at  y = 0 (plane  of  symmetry).  The  following  is  a brief  description 
of  the  main  program  and  subroutines  (listed  alphabetically)  shown  in 
Figure  7. 

EPIC -3  This  is  the  main  program  which  calls  subroutines  MATE, 

GEOM,  START  and  LOOP.  For  restart  calculation,  the 
program  goes  directly  to  subroutine  LOOP  (56  cards). 

BREAK  This  subroutine  checks  for  fracture  (20  cards). 


BRICK  This  subroutine  generates  a series  of  composite  brick 

elements,  each  composed  of  six  tetrahedron  elements 
(A,  B,  C,  D,  E,  F)  as  shown  in  Figure  8 (49  cards). 

DETERM  This  subroutine  calculates  the  determinate  of  an  array 
for  volume  calculations  (50  cards). 
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ELEG 

E SHAPE 

GEOM 

GPLOT 

LOAD 

LOOP 

MATL 

NODEG 

NSHAPE 

RECALL 


This  subroutine  generates  data  for  a single  tetrahedron  element 
and  places  this  initial  data  into  the  appropriate  internal  file 
(43  cards). 

This  subroutine  generates  elements  for  various  special  shapes 
composed  of  rod  geometry,  rounded-nose  geometry,  conical- 
nose  geometry  and  flat-plate  geometry  (1269  cards). 

This  subroutine  reads,  generates  and  prints  the  initial 
geometry  (451  cards). 

This  subroutine  writes  data  on  TapeS  for  geometry  plots 
(52  cards). 

This  subroutine  calculates  axial  forces,  shear  forces  and 
bending  moments  in  a slender  projectile  (69  cards). 

This  subroutine  calculates  the  equations  of  motion,  strains, 
strain  rates,  stresses  and  concentrated  nodal  forces.  It 
also  prints  selected  node  and  element  data  (520  cards). 

This  subroutine  reads,  calculates  and  prints  the  material 
properties  (156  cards). 

This  subroutine  generates  a row  of  nodes  (69  cards). 

This  subroutine  generates  nodal  geometry  for  various 
special  shapes  composed  of  rod  geometry,  rounded-nose 
geometry,  conical-nose  geometry  and  flat-plate  geometry 
(581  cards). 

This  subroutine  reads  data  from  Tape2  for  a restart  run 
(88  cards). 
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SAVE 


SDATA 


SLIDE 


START 


STRESS 


TETRA 


WEDGEl 


WEDGE2 


This  subroutine  writes  data  on  Tape2  for  a restart  run 
(76  cards)* 

This  subroutine  calculates,  prints  and  saves  (Tape4)  system 
data  such  as  eg  positions,  kinetic  energies,  and  linear  and 
angular  momenta  for  the  projectile,  target  and  total  system 
(229  cards). 

This  subroutine  calculates  the  equations  of  motion  of  the  slave 
nodes.  If  there  is  interference,  ii  places  the  slave  node  on 
the  master  surface  and  transfers  the  lost  momenta  of  the 
slave  nodes  to  the  appropriate  master  nodes  (2  56  cards). 

This  subroutine  reads  the  initial  dynamic  conditions  and  sets 
the  variables  equal  to  the  conditions  which  exist  at  impact 
(108  cards). 

This  subroutine  calculates  element  stresses  which  consist  of 
elastic  stresses,  plastic  and  viscous  deviator  stresses, 
hydrostatic  pressure  and  artificial  viscosity  (89  cards). 

This  subroutine  generates  a series  of  tetrahedron  elements 
(33  cards). 

This  subroutine  generates  a series  of  composite  wedge  ele- 
ments, each  composed  of  three  tetrahedron  elements  contain- 
ing nodes  Nl,  N3,  N4,  N5,  N7,  N8,  (Elements  A,  B,  C),  as 
shown  in  Figure  8 (40  cards). 

This  subroutine  generates  a series  of  composite  wedge  ele- 
ments, each  composed  of  three  tetrahedron  elements  contain- 
ing nodes  N 1,  N2,  N3,  N5,  N6,  N7,  (Elements  D,  E,  F),  as 
shown  in  Figure  8 (40  cards). 


The  plotting  program  reads  from  Tape3  (written  by  subroutine  GPLOT)  and 
generates  a two-dimensional,  cross-sectional  plot  of  the  x-z  plane  at  y = 0. 


3,  2 PROGRAM  LOGIC  AND  STORAGE  ALLOCATIONS 

Program  logic  and  storage  allocation  are  shown  schematically  in  Figure  9. 
The  existing  EPIC-3  code  is  dimensioned  for  a maximum  of  2000  nodes,  as 
shown.  The  element  data  are  stored  on  two  disk  files  in  blocks  of  200  ele- 
ments each.  Since  only  one  block  of  element  data  is  in  core  at  any  specific 
time,  there  is  no  limit  to  the  number  of  elements.  Relatively  small  storage 
allocations  are  made  for  material  properties  and  sliding-surface  data  which 
identify  the  master  nodes  and  slave  nodes.  The  EPIC~3  code,  with  a 
capacity  of  2000  nodes,  requires  69K  storage  in  a Honeywell  6080  computer. 
There  are  generally  about  four  to  five  times  as  many  elements  as  nodes, 
so  this  version  of  the  program  could  provide  solutions  for  problems  involv- 
ing up  to  2000  nodes  and  about  10,  000  elements.  It  should  also  be  noted  that 
many  computers  (including  the  Honeynvell  6080)  have  storage  capacities 
significantly  greater  than  69K.  Therefore,  by  increasing  the  dimensions  of 
the  node  arrays,  solutions  to  much  larger  problems  can  be  obtained  with  the 
existing  code. 

Most  of  the  computations  occur  in  two  large  loops,  the  node  loop  and  the 
element  loop.  The  primary  activity  for  the  node  loop  is  to  apply  the  equa- 
tions of  motion  to  all  nodes  using  the  updated  accelerations  from  the  con- 
centrated forces  which  were  accumulated  during  the  previous  element  loop. 

The  element  loop  operates  on  each  element  in  a consecutive  order  beginning 
with  the  first  element  in  the  first  block,  and  ending  with  the  last  element 
in  the  last  block.  The  important  activities  in  the  element  loop  are: 
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P'igure  9.  EPIC -3  j-rogram  Logic  and  Storage 


• Bring  block  of  element  data  from  Pale  A when  computations  on 
previous  block  are  complete, 

• Obtain  updated  element  data  such  as  strains  and  strain  rates. 
Since  the  strains  and  strain  rates  are  dependent  on  the  current 
displacements  and  velocities  of  the  nodes,  pointers  are  required 
to  obtain  these  data  from  the  node  arrays  \vhich  are  always  in 
the  central  memory  core. 

• Obtain  stresses  in  the  elements  using  the  updated  strains  and 
strain  rates.  This  activity  requires  a pomter  to  the  material 
arrays. 
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• Calculate  the  concentrated  nodal  forces  from  the  updated 
stresses,  and  distribute  them  to  the  appropriate  node  arrays 
with  the  same  pointers  used  to  obtain  the  displacements  and 
velocities. 

• When  the  entire  block  of  element  data  has  been  processed, 
write  the  updated  data  on  File  B and  bring  the  next  block  of 
element  data  into  core. 

It  should  be  noted  that  all  blocks  (and  therefore  all  elements)  are  processed 
before  the  next  node  loop  begins.  Also,  for  the  element  loop  computations 
during  the  next  integration  time  increment,  disk  Files  A and  B are  inter- 
changed. 


3.  3 NODE  AND  ELEMENT  ARRAYS 

There  are  16  node  variables  which  are  contained  in  arrays  and  23  element 
variables  contained  in  arrays  as  shown  in  Figure  9.  The  node  arrays  are 
described  as  follows: 

NODE  Identifies  the  node  number  such  that  it  is  not  necessary  to 

have  a consecutive  node  numbering  system 

AMASS  The  total  mass  of  the  node 

PMASS  The  mass  of  the  node  contributed  by  the  projectile 

XI  The  initial  x coordinate  of  the  node 

YI  The  initial  y coordinate  of  the  node 

ZI  The  initial  z coordinate  of  the  node 
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X 


The  displaced  x coordinate  of  the  node 


Y The  displaced  y coordinate  of  the  node 

Z The  displaced  z coordinate  of  the  node 

XDOT  The  x velocity  of  the  node 

YDOT  The  y velocity  of  the  node 

ZDOT  The  z velocity  of  the  node 

PX  The  X direction  concentrated  force  on  the  node 

PY  The  y direction  concentrated  force  on  the  node 

PZ  The  z direction  concentrated  force  on  the  node 

IFIX  Identifies  X,  Y,  Z restraints  on  the  node.  It  also  identifies 

a slave  node  and  indicates  if  it  is  free  or  sliding  on  the  master 
surface. 

The  element  arrays  are  described  as  follows: 

NEL  Identifies  the  element  number 

NODE  1 Node  number  of  node  i (see  Figure  2) 

NODE  2 Node  number  of  node  j 

NODE  3 Node  number  of  node  m 
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NODE  4 


MAT 

BI 

BJ 

BM 

BP 

Cl 

CJ  — 

CM 

CP 

DI 

DJ 

DM 

DP__ 

VOL 

DVOL 

ES 

EW 

ICHECK 


Node  number  of  node  p 
Material  number 


Initial  geometry-dependent  constants  defined  in  Equations 
(6c)  through  (6e). 


Initial  volume 

Volumetric  strain  defined  in  Equation  (13) 

Total  internal  energy  per  original  unit  volume  defined  in 
Equation  (56) 

Heat  energy  per  original  unit  volume 

Identifies  state  of  material  (elastic,  plastic,  fractured) 
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SECTION  IV 

PROGRAM  USER  INSTRUCTIONS 


The  required  EPIC-3  input  data  for  various  conditions  are  summarized  in 
Figure  10  and  described  in  subsections  4.  1 through  4.4.  Output  data  are 
described  in  subsection  4.  5.  Succeeding  subsections  deal  with  diagnostics, 
estimated  central  processor  time  requirements,  and  central  memory 
storage  requirements  and  alterations. 


4.  1 INPUT  DATA  FOR  AN  INITIAL  RUN 


An  initial  run  is  one  which  requires  input  data  for  the  initial  geometry  and 
impact  conditions  at  time  = 0.  The  descriptions  which  follow  are  for  the 
input  data  in  Figure  10.  Consistent  units  must  be  used. 


• Identification  Card  (215,  FIO.  0)  - 

CASE  = Case  number  for  run  identification 
CYCLE  ==  0 for  initial  run 

CPMAX  = Central  Processor  time  (hours)  at  which  the  results 
will  be  written  into  the  restart  tape  and  the  run  will 
stop.  This  feature  will  be  b3rpassed  if  CPMAX  = 0. 


• 24  Material  Cards  (5E15.8)  - Each  of  the  24  material  cards  provides 

a specific  material  characteristic  for  five  materials.  For  instance, 
the  first  card  describes  the  density  of  the  materials.  The  density  of 
material  1 is  entered  in  columns  1-15,  the  density  of  material  2 in 
columns  16-30,  etc.  The  second  card  contains  the  specific  heats  for 
the  five  materials,  etc.  It  is  only  necessary  to  describe  the  mater- 
ials used  for  the  analysis.  If  materials  1 and  3 are  used  for  a specific 


INITIAL  RUN  INPUT  DATA 


IDENTIFICATION  CARD  (215,  FIO.O) 


24  MATERIAL  CARDS  (5E15.8) 


MATLl  1 MATL2 

MATL3  MATL4 

MATL5 

PROJECTILE  SCALE  / SHIFT/ ROTATE  CARD  (6F10.0) 

X SCALE  YSCALE  Z SCALE 

XSHIFT  1 ZSHIFT  ROTATE 

PROJECTILE  NODE  DATA  CARDS  - AS  REQUIRED  (315,  2X,  311,  7F8.0) 


NN 

INC 

T 

XI 



Y1 

Z1 

XN 

YN 

ZN 

EXPAND 

IX, lY,  IZ 


BLANK  CARD 


PROJECTILE  NODE  SPECIAL  SHAPES  - AS  REQUIRED  (DESCRIPTION  FOLLOWS) 
BLANK  CARD  I 


TARGET  SCALE/ SHIFT/ ROTATE  CARD  (6F10.0) 


XSCALE 

YSCALE 

Z SCALE 

XSHIFT 

ZSHIFT 

ROTATE 

TARGET  NODE  DATA  CARDS  - AS  REQUIRED  (315,  2X,  311,  7F8.0) 


N1 

NN 

INC 

i; 

n 

XI 

Y1 

Z1 

XN 

YN 

ZN 

EXPAND 

IX, lY,  IZ 


BLANK  CARD  \ 

TARGET  NODE  SPECIAL  SHAPES  - AS  REQUIRED  (DESCRIPTION  FOLLOWS) 
BLANK  CARD  j 

PROJECTILE  ELEMENT  DATA  CARDS  - AS  REQUIRED  (1315) 


ELEN 

INC 

MATL 

N1 

N2 

N3 

N4 

N5 

N6 

N7 

N8 

l«tli]< 

WSm 

BLANK  CARD 


PROJECTILE  ELEMENT  SPECIAL  SHAPES  - AS  REQUIRED  (DESCRIPTION  FOLLOWS) 
BLANK  CARD  I 


Figure  10.  EPIC -3  Input  Data  Summary 
(1  of  3) 
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TARGET  ELEMENT  DATA  CARDS  “ AS  REQUIRED  (1315) 

I ELEl  I ELEN  I fjljc^WTlj  N1  [ N2  | N3  f N^  N5  [ No  7 


T“‘^:riNOi5ri 

N/  No  i [Mr 


BLANK  CARD 


TARGET  ELEMENT  SPECIAL  SHAPES  - AS  REQUIRED  (DESCRIPTION  FOLLOWS) 
BLANK  CARD  ] 


SLIDE /LOAD  CARD  (515) 
InSLIdI  I R IG  NTOpTnBOT  [nR  1 


ixjuw  iivio  wiur  nDui  i nohw ^ v- ' 

i ^ 1 1 ! 

SLIDING  SURFACE  IDENTIFICATION  CARDS  - AS  REQUIRED  (415) 


TYPE  NMX  NMY  NSLA#?^^ 


A4ASTER  NODES  - AS  REQUIRED  (1615)  NMY  ROWS  WITH  NMX  NODES  PER  ROW 

START  EACH  ROW  WITH  NEW  CARD 

iMi  1 1M2  i '1  ^ I I ^ I I 1 ^ r 


IM15  IM16 


SLAVE  NODES  - AS  REQUIRED  (1615) 

] I 1 I ■ r 

isi  Ms2  i ! i 


INITIAL  VELOCITY  CARD  (6F10.0) 


) IS16 


PXDOT  PYDOT  r PZDOT  T TXDOT  ^ TYDOT  ! TZDo"t 

L ^ 1 

INTEGRATION  TIME  INCREMENT  CARD  (5rl0,0) 


DTI  DTMAX  DTMiN 


TMAX  f 


DATA  OUTPUT  CARDS  - AS  REQUIRED  (2F10.0,  415) 

T 1 ME  I ECHECK  TTs  ^^LOAdH  SA  VENP  LOtKP 


RESTART  RUN  INPUT  DATA 


IDENTIFICATION  CARD  (215,  FIO.O) 

[case  cycle  CP  max 

INTEGRATION  TIME  INCRE/^tNT  CARD  (3F10,0) 


DTMIN  STF 


DATA  OUPUT  CARDS  - AS  REQUIRED  (2F10.0,  415) 
TIME  r ECHECK  flSYsTlLOApi  ISAVeIipLOT^ 


Figure  10.  EPIC'S  Input;  Ijata  Suminarv 
(2  ol  3) 
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SPECIAL  SHAPES  INPUT  DATA 


3 NODE  CARDS  FOR  EACH  ROD  SHAPE  (415,  3F10. 0 / 5F10. 0 / 5F10. 0) 


1 

N1 

NRING 

NPLN 

1 ZTOP 

ZBOT 

RTOPl 

RTOP3 

RTOP4 

RBOTl 

R60T2 

RBOT3  1 RBOT4  | RBOT5 

3 NODE  CARDS  FOR  EACH  NOSE  SHAPE  (315,  5X,  FIO.O  / 5F10.0  / 5F10.0) 
INOSE  - 2--»-ROUNDED  NOSE,  INOSE  - 3 -^CONICAL  NOSE 


INOSE 

RTOPl 

RTOP2 

RTOP4 

RTOP5 

ZMINl 

ZM1N3 

ZMIN4 

3 NODE  CARDS  FOR  EACH  FLAT  PLATE  SHAPE  (215  / 615,  3F10.0  / 6F10.0) 

4 

N1 

i 

1 

1 

1 

p 

1 

1 

NX 

NY 

NZ  1 IDY  I IDZ  lY 

1 

X-EXPAND 

XI 

Y1 

Z1 

XN 

YN 

ZN 

ELEMENT  CARD  FOR  EACH  ROD  SHAPE  (1015) 


1 

N1 

NRING 

NLAY 

ELEl 

Ml 

M2 

M3 

M4 

M5 

ELEMENT  CARD  FOR  EACH  NOSE  SHAPE  (315,  5X,  615) 

INOSE  - 2— ROUNDED  NOSE,  INOSE  - 3— --CONICAL  NOSE 


INOSE  N1 


NRING 


EL£1 


Ml 


M2 


M3 


IVM 


ELEMENT  CARD  FOR  EACH  FLAT  PLATE  SHAPE  (315,  5X,  515) 


Figure  10.  EPIC -3  Input  Data  Summary 
(3  of  3) 
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analysis,  all  the  material  1 data  must  be  specified  in  columns  1-15 
and  all  the  material  3 data  in  columns  31-45,  Columns  16-30,  46-60, 
and  61-75,  representing  materials  2,  4 and  5,  can  be  left  blank. 


Card  1 
Card  2 
Card  3 
Card  4 
Card  5 
Card  6 
Card  7 
Card  8 


Card  9 
Card  10 
Card  11 
Card  12 
Card  13 
Card  14 
Card  15 

Card  16 


Density  (mass /volume) 

Specific  heat  (work/mass /degree) 

Modulus  of  elasticity  (force /area) 

Poisson’s  ratio 

Absolute  viscosity  (force-time/area) 

Yield  stress  (force/area) 

Ultimate  stress  (force/area) 

Strain  at  which  the  ultimate  stress  is  achieved.  Should  be 
consistent  with  the  strain  definition  in  Equation  (14). 

Stress  varies  linearly  between  the  strain  at  the  yield 
stress  and  the  strain  at  the  ultimate  stress. 

The  stress  for  larger  strains  is  equal  to  the  ultimate 
stress  until  fracture  occurs. 

Maximum  negative  hydrostatic  pressures  (force /area) 

Strain  rate  effect  constant,  Ch , in  Plquation  (31) 

Pressure  effect  constant,  C^,  in  Equation  (31) 

Pressure  effect  constant,  C^,  in  Equation  (31) 

Temperature  effect  constant,  C.,,  in  Equation  (31) 

Temperature  effect  constant,  in  Equation  (31) 

Hydrostatic  pressure  constant,  , in  Equation  (32) 
(force/area) 

Hydrostatic  pressure  constant,  K.^,  in  Equation  (3  2) 
(force/  area) 


Card  17 

Card  18 
Card  19 
Card  20 
Card  21 
Card  22 
Card  23 


Card  24 


Hydrostatic  pressure  constant,  K^,  in  Equation  (32) 

(force /area) 

Gruheisen  coefficient,  r a Eqtiation  (32) 

Linear  artificial  viscosity  coefficient,  C^,  in  Equation  (33) 

2 . 

Quadratic  artifical  viscosity  coefficient,  in  Equation  (33) 


Equivalent  strain 
Volumetric  strain 


if  both  are  exceeded,  the  element  is  totally 
failed  and  produces  no  stresses  or  pressures. 


Equivalent  strain  for  shear  and  tensile  failure  only.  Positive 
hydrostatic  pressure  and  viscosity  stress  capability  remain. 
If  negative,  the  material  behaves  like  a liquid. 


Initial  temperature  (degrees) 


• Projectile  Scale / Shift/Rotate  Card  (6F1Q.  0)  - 

XSCALE  = Factor  by  which  the  x coordinates  of  all  projectile 
nodes  are  multiplied.  Applied  after  the  coordinate 
shifts  (XSHIFT,  ZSHIFT)  described  later. 


YSCALE  “ Factor  by  which  the  y coordinates  are  multiplied. 


ZSCALE  = Factor  by  which  the  z coordinates  are  multiplied. 


XSHIFT  ” Increment  added  to  the  x coordinates  of  all  projectile 
nodes  (length).  Applied  before  the  scale  factors 
(XSCALE,  YSCALE,  ZSCALE). 

ZSHIFT  = Increment  added  to  the  z coordinates  (length). 


ROTATE  “ Rotation  about  the  y axis  (at  x=z  = 0)  of  all  projectile 
nodes  (degrees).  Applied  after  the  coordinate  shifts 
(XSHIFT,  ZSHIFT),  and  the  scale  factors  (XSCALE, 
YSCALE,  ZSCALE), 
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• Projectile  Node  Data  Cards  (315,  2X,  311,  7F8,  0)  - The  projectile  node 
data  cards  are  provided  as  required  for  the  projectile  nodes.  If  a 
node  is  at  the  interface  of  the  projectile  and  the  target,  and  contains 
mass  from  both  the  projectile  and  the  target,  it  must  be  included  with 
the  projectile  nodes.  The  node  numbers  (Nl*  - * NN)  must  not  exceed 
the  dimension  of  the  node  arrays,  and  they  need  not  be  numbered 
consecutively  or  in  increasing  order.  End  with  a blank  card. 


N1 

First  node  in  a row  of  nodes. 

NN 

Last  node  in  a row  of  nodes.  Leave  blank  if  a single 
node  (not  a row)  is  entered. 

INC 

Node  number  increment  between  adjacent  nodes.  If 
left  blank,  program  will  set  INC  - 1. 

IX 

1 will  restrain  all  nodes  (Nl-  • ^ NN)  in  the  x direction. 

No  restraint  if  left  blank. 

lY 

1 will  restrain  nodes  in  y direction. 

IZ 

=: 

1 will  restrain  nodes  in  z direction 

XI 

X coordinate  of  node  N1  (length) 

Y1 

y coordinate  of  node  N1  (length) 

Z1 

z coordinate  of  node  N1  (length) 

XN 

= 

X coordinate  of  node  NN  (length).  Leave  blank  if  a 
single  node  is  entered. 

YN 

rr 

y coordinate  of  node  NN  (length) 

ZN 

z coordinate  of  node  NN  (length) 

EXPAND 

Factor  by  which  the  distance  between  adjacent  nodes 
changes.  For  example,  if  the  distance  between  the 

first  two  nodes  is  A^,  the  distance  between  the  second 
and  third  nodes  is  a 9 ■“  A-j^  ' EXPAND.  A summary  of 

49 


relative,  spacing  between  the  first  two  nodes  and  the 
last  two  nodes  is  given  in  Figure  11.  If  left  blank, 
the  program  will  set  EXPAND  = 1.  0 for  uniform 
spacing. 

• Projectile  Node  Special  Shapes  - These  are  described  later.  End  with 
a blank  card. 

• Target  Scale /Shift /Rotate  Card  (6F1Q.  0)  - Same  as  Projectile  Scale/ 
Shift /Rotate  Card  except  it  applies  to  the  target  nodes. 

• Target  Node  Data  Cards  (315.  2X,  311,  78F.  0)  - Same  as  Projectile  Node 
Data  Cards  except  these  apply  to  the  target  nodes.  If  there  are  inter- 
face nodes  which  contain  mass  from  both  the  projectile  and  the  target, 
these  nodes  are  included  with  the  projectile  nodes  and  must  not  be 
included  here.  End  with  a blank  card. 

• Target  Node  Special  Shapes  - These  are  described  later.  End  with  a 
blank  card. 

• Projectile  Element  Data  Cards  (1315)  - The  Projectile  Element  Data 
Cards  are  provided  as  required  for  the  projectile  elements.  For  this 
immediate  discussion  it  will  be  assumed  that  the  elements  are  entered 
as  a series  of  composite  brick  elements,  each  containing  six  individual 
elements  as  shown  in  Figure  8.  Following  this  immediate  discussion 
will  be  an  example  and  instructions  for  generating  a series  of 
individual  elements.  End  with  a blank  card, 

ELEl  First  individual  element  in  the  first  composite  brick 

element  in  a series  of  composite  brick  elements. 

The  six  elements  in  a composite  brick  element  are 
numbered  consecutively. 
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NUMBER 


ELEN 


First  individual  element  in  the  last  composite  brick 
element  in  a series  of  composite  brick  elements. 

Note  that  this  is  the  first,  rather  than  last,  individual 
element  in  the  last  composite  brick  element.  Leave 
blank  if  a single  composite  element  is  entered. 

ELE  INC  = Element  number  increment  between  the  first  individual 
element  of  successive  composite  brick  elements. 

Since  each  composite  brick  element  contains  six  indivi- 
dual elements,  this  increment  must  not  be  less  than  six. 
If  left  blank  the  program  will  set  ELE  INC  = 6. 

MATL  = Material  number  (1,  2,  3,  4,  or  5)  of  the  elements.  If 
left  blank,  the  material  number  from  the  previous  ele- 
ment data  card  will  be  used. 

N1  - N8  = Node  numbers  of  the  first  composite  brick  element  as 
shown  in  Figure  8.  Nodes  Nl,  N2,  N3  and  N4,  and 
nodes  N5,  N6,  N7  and  N8,  are  counterclockwise  when 
looking  from  Nl  to  N5. 

NODE  INC  = The  node  number  increment  added  to  the  node  numbers 
of  the  previous  composite  brick  element  for  the  next 
composite  brick  element.  If  left  blank,  the  program 
will  set  NODE  INC  = 1. 

An  example  of  input  data  for  composite  brick  elements  is  shown  in 
Figure  12,  In  the  upper  left  corner  it  can  be  seen  that  there  are  four 
rows  of  nodes  (12-18,  22-28,  32-38,  42-48)  which  are  arranged  to 
contain  three  composite  brick  elements.  If  ELEl  * 101,  ELEN  * 121 
and  ELE  INC  = 10,  then  the  first  composite  brick  contains  elements 
101-106,  the  second  contains  111-116,  and  the  third  contains  121-126, 
The  first  composite  brick  is  defined  by  nodes  Nl  = 12,  N2  = 22,  N3  - 32, 
N4  = 42,  N5  = 14,  N6  = 24,  N7  = 34  and  N8  = 44,  Note  that  N1-N4  and 
N5-N8  are  counterclockwise  when  looking  from  Nl  to  N5,  The  six 
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individual  elements  are  generated  according  to  the  arrangement  and 
order  (A,  B,  C,  D,  E,  F)  shown  in  Figure  8.  The  node  numbers  for 
each  successive  brick  are  simply  NODE  INC  = 2 greater  than  those  of 
previous  brick.  For  the  second  brick  for  instance,  N1  = 12  + 2 = 14, 

N2  = 22  + 2 = 24,  N3  = 32  + 2 = 34,  N4  = 42  + 2 = 44,  N5  = 14  + 2 = 16, 

N6  = 24  + 2 = 26,  N7  = 34  + 2 = 36  and  N8  = 44  + 2 = 46. 

It  is  possible  to  generate  a series  of  individual  tetrahedron  elements  by 
letting  N1-N4  be  the  nodes  of  the  first  element,  where  N 1,  N2,  and  N3 
are  counterclockwise  when  viewed  from  N4,  This  option  is  exercised 
when  N5-N8  are  left  blank.  For  this  option,  if  ELE  INC  is  left  blank, 
the  program  will  set  ELE  INC  = 1,  It  is  also  possible  to  generate  a 
series  of  composite  wedge  elements,  each  containing  three  individual 
tetrahedron  elements.  The  three  elements  in  a composite  wedge  ele- 
ment are  numbered  consecutively.  If  N2  and  N6  are  left  blank,  the  first 

three  tetrahedron  elements  (A,  B,  C)  are  defined  by  nodes  N 1,  N3,  N4, 

N5,  N7  and  N8  as  shown  in  Figure  8.  Likewise,  if  N4  and  N8  are  left 
blank,  the  first  three  elements  (D,  E,  F)  are  defined  by  nodes  Nl,  N2, 
N3,  N5,  N6  and  N7,  For  the  composite  wedge  element  options,  if 
ELE  INC  is  left  blank,  the  program  sets  ELE  INC  = 3,  which  is  the 
minimum  acceptable  increment. 

• Projectile  Element  Special  Shapes  - These  are  defined  later.  End 
with  a blank  card. 

• Target  Element  Data  Cards  (1315)  - Same  as  Projectile  Element  Data 
Cards.  End  with  a blank  card. 

• Target  Element  Special  Shapes  - These  are  defined  later.  End  with 
a blank  card, 

• Slide /Load  Card  (515)  - 

NS  LID  = The  number  of  sliding  surfaces.  Currently  limited  to  a 

maximum  of  two.  This  does  not  include  the  rigid  sliding 
surface  described  next. 
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IRIG 


1 gives  a rigid  frictionless  surface  on  the  positive 
side  of  the  plane  described  by  z = 0.  If  the  equations 
of  motion  cause  a node  to  have  a negative  z coordinate, 
the  node  coordinate  is  redefined  and  set  at  z = 0.  If 
IRIG  = 0 this  option  is  not  used. 

The  next  three  variables  (NTOP,  NBOT,  NRING)  are  used  to  identify 
the  nodal  geometry  for  internal  load  calculations  in  a slender  projectile. 
If  this  option  is  used,  the  nodal  geometry  must  be  consistent  with  the 
arrangement  of  nodes  generated  by  the  special  shapes  generator  for  rod 
geometry.  Leave  blank  if  it  will  not  be  used. 

NTOP  = The  node  number  of  the  centerline  node  at  the  top  free 

end  of  the  cylindrical  projectile. 

NBOT  = The  node  number  of  the  centerline  node  at  the  lower 

end  of  the  cylindrical  portion  of  the  projectile, 

NRING  = The  number  of  rings  of  elements  in  the  cylindrical 

projectile.  This  will  be  described  in  the  discussion  of 
special  shapes. 

• Sliding  Surface  Identification  Cards  (415)  - Each  sliding  surface  con- 
tains one  identification  card  and  cards  (as  required)  describing  the 
master  nodes  and  the  slave  nodes.  If  there  are  no  sliding  surfaces 
(NSLID  = 0),  the  Initial  Velocity  Card  (described  later)  follows  the  pre- 
viously described  Slide/ Load  Card.  If  there  are  two  sliding  surfaces, 
the  Sliding  Surface  Identification  Card  for  the  second  sliding  surface 
follows  the  slave  node  data  for  the  first  sliding  surface, 

TYPE  = 1 gives  the  slave  surface  above  the  master  surface 

relative  to  the  z axis,  TYPE  = 2 gives  slave  surface 
below  the  master  surface  relative  to  the  z axis.  In 
both  cases  the  z coordinates  of  the  master  surface 
must  be  a single -valued  function  of  x and  y as  dis- 
cussed in  subsection  2.  6. 
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• Data  Output  Cards  (2F1Q,  0,  415)  - These  cards  are  used  to  specify 

various  forms  of  output  data  at  selected  times,  and  the  last  card  must 
be  for  a time  greater  than  TMAX,  even  though  output  will  not  be  pro- 
vided for  that  specific  time. 

TIME  = Time  at  which  output  data  will  be  provided. 

ECHECK  = A code  which  governs  the  printed  output.  The  three 

following  options  are  provided: 

1)  If  ECHECK  is  greater  than  1000.  , the  individual 
node  data  and  element  data  will  not  be  printed. 

Only  system  data  such  as  eg  positions,  momenta, 
kinetic  energies  and  average  velocities  are  pro- 
vided for  the  projectile,  target  and  entire  system. 

2)  If  ECHECK  is  positive  (includes  0)  and  less  than 
1000. , the  system  data  and  individual  node  data 
will  be  printed.  Individual  element  data  will  be 
printed  for  all  elements  which  have  an  equivalent 
strain  [Equation  (14)]  equal  to  or  greater  than 
ECHECK.  For  example,  if  ECHECK  = 0. , all 
element  data  will  be  printed.  If  ECHECK  = 0.  5, 
only  those  elements  with  equivalent  strains  equal 
to  or  greater  than  0.  5 will  have  data  printed.  If 
ECHECK  = 999.  , no  element  data  will  be  printed. 

3)  If  ECHECK  is  negative,  the  system  data  and 
individual  node  data  will  be  printed.  The  individual 
element  data  will  be  printed  only  for  those  ele- 
ments which  are  plastic  or  failed.  Data  for  the 
elastic  elements  will  not  be  printed. 
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ISYS 


1 will  write  system  data  on  Tape4.  (Notei  A program 
has  not  yet  been  written  to  read  and  plot  these  results. ) 


ILOAD  = 1 will  calculate  internal  loads  in  a cylindrical  projec- 

tile. This  option  can  only  be  used  if  NTOP,  NBOT 
and  NPING  are  properly  defined  in  the  Slide /Load  Card. 

ISAVE  = 1 will  write  results  on  Tape2  for  a possible  restart  run. 

IPLOT  = 1 will  write  results  in  TapeS  for  a possible  plot. 

With  the  exception  of  the  special  shapes,  this  completes  the  description 
of  the  input  data  for  an  initial  run.  Four  types  of  special  shapes  can 
be  generated  for  either  the  projectile  or  the  target.  The  node  and  ele- 
ment data  for  these  special  shapes  are  entered  individually  in  the  loca- 
tions identified  in  Figure  10.  The  node  data  for  each  special  shape  re- 
quire three  cards,  and  the  element  data  require  one  card  for  each 
special  shape.  A description  of  the  special  shapes  data  follows, 

• 3 Node  Cards  for  Each  Rod  Shape  (415.  3F10.  0/5F10.  0/5F10.  0)  - The 

rod  geometry  descriptions  are  given  in  Figure  13,  The  rod  is  always 
generated  in  a vertical  position  about  the  z axis.  Only  half  the  rod  is 
generated  as  shown,  and  restraints  are  provided  normal  to  the  plane  of 
symmetry.  The  rotation  of  the  rod  for  oblique  impact  is  obtained  with 
a Scale/Shift/Rotate  Card. 

1 = Identification  number  for  rod  geometry. 

N1  = The  node  number  of  the  centerline  node  on  the  top  end 

of  the  rod.  The  nodes  of  the  rod  are  numbered  con- 
secutively beginning  with  Nl.  It  can  be  seen  that  the 
first  ring  on  the  top  plane  is  defined  by  nodes  N2-N6, 
in  a counterclockwise  direction.  The  second  ring  is 
defined  by  nodes  N7-N15,  the  third  ring  by  nodes 
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NUMBER 
OF  RINGS 
(NRING) 

NODES 
PER  PLANE 
(NPP) 

ELEMENTS 
PER  LAYER 
(EPL) 

1 

6 

12 

2 

15 

48 

3 

24 

96 

4 

37 

156 

5 

50 

228 

N16-N24,  the  fourth  ring  by  nodes  N25-N37,  and  the 
fifth  ring  by  nodes  N38-N50,  The  node  numbers  in  the 
plane  of  nodes  below  the  top  plane  are  numbered  in  the 
same  order  as  the  top  plane  except  each  number  is  NPP 
(Nodes  Per  Plane)  greater  than  the  node  above  it. 

Values  for  NPP  are  dependent  on  the  number  of  rings 
and  are  shown  in  Figure  11.  (Example:  If  N1  = 101 
and  there  are  four  rings,  the  centerline  node  in  the 
second  plane  is  N1  + NPP  = 101  + 37  = 138. ) With  this 
technique  it  is  possible  to  identify  readily  the  appro- 
priate surface  nodes  for  slave  node  designation  and  the 
bottom  centerline  node,  NBOT,  for  the  internal  loads 
calculat  ions. 

NRING  = The  number  of  rings  (1,  2,  3,  4 or  5)  in  the  rod. 

NPLN  = The  number  of  cross-sectional  planes  of  nodes  in  the 

rod. 

ZTOP  = The  z coordinate  of  the  top  of  the  rod. 

ZBOT  = The  z coordinate  at  the  bottom  of  the  rod, 

EXPAND  = Factor  by  which  the  vertical  distance  between  adjacent 
node  changes.  Factor  applies  from  top  to  bottom. 

Same  as  described  for  the  Projectile  Node  Data  Cards. 

RTOPl  - RTOP5  = Radii  of  the  rings  at  the  top  of  the  rod, 

RBOTl  - RBOT5  Radii  of  the  rings  at  the  bottom  of  the  rod. 

Note:  F it  is  not  possible  to  describe  the  node  geometry  of  the 

rod  with  a single  shape,  it  is  possible  to  use  multiple 
shapes  to  form  a single  rod.  The  nodes  must  be  numbered 
consecutively,  and  NRING  must  be  the  same  for  all  the 
individual  rod  shapes. 
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3 Node  Cards  for  Each  Nose  Shape  (315,  5X.  FIO.  0/5F10.  0/5F10,  0)  - 
The  nose  geometry  descriptions  are  given  in  Figure  14*  The  nose 
shape  is  always  generated  in  a vertical  position  (pointed  downward) 
about  the  z axis.  Only  half  the  nose  shape  is  generated  as  shown,  and 
restraints  are  provided  normal  to  the  plane  of  symmetry.  The  node 
geometry  for  the  flat  plane  at  the  rod  interface  is  not  generated  with  the 
nose  generator  and  must  therefore  be  generated  with  the  rod  generator. 
The  number  of  rings  must  be  identical  for  the  rod  and  the  nose.  The 
surface  nodes  are  also  identified  in  Figure  14  so  they  can  be  designated 
as  slave  nodes  if  necessary. 


INOSE 


2 identifies  a rounded  nose  shape.  If  the  length  of  the 
nose  is  equal  to  the  radius,  a spherical  shape  is  gener- 
ated. When  the  length  is  not  equal  to  the  radius,  the 
axial  coordinates  are  scaled,  and  the  various  radii  are 
not  changed. 


INOSE 


3 identifies  a conical  nose  shape. 


The  node  number  of  the  centerline  node  at  the  top  of  the 
nose  at  the  rod  interface.  This  is  identical  to  the  bottom 
centerline  node  of  the  rod.  This  node  is  not  the  same  as 
node  N1  used  for  the  rod  shape.  This  node  (and  the  other 
rod  interface  nodes)  is  not  included  in  the  nose  generator 
and  must  be  generated  with  the  rod  generator.  It  is  used 
for  reference  such  that  the  nodes  can  number  in  a consecu- 
tive manner. 


NRING  = Number  of  rings  in  the  nose  and  the  rod. 


ZTOP 


The  z coordinate  at  the  top  of  the  nose.  This  is  identical 
to  ZBOT  for  the  rod  shape  at  the  rod  interface. 
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ROUNDED  NOSE  (INOSE  = 2) 


CONICAL  NOSE  (INOSE  = 3) 


* DOES  NOT  INCLUDE  NODES  AT  ROD  INTERFACE 


Figure  14.  Geometry  of  Special  Nose  Shapes 
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RTOPl  - RTOP5  = Radii  of  the  rings  at  the  top  of  the  nose  shape. 

These  are  generally  equal  to  RBOTl  - RBOT5 
for  the  rod  shape  at  the  rod  interface. 


ZMINl  - ZMIN5  = Minimum  coordinates  for  the  rings  in  the  nose 

shape  (i.  e.  , z coordinates  of  the  nodes  which 
lie  on  the  z axis). 


plate  descriptions  are  given  in  Figure  15.  In  all  cases,  the  lines  con- 
necting the  adjacent  corner  nodes  are  parallel  to  one  of  the  three  pri- 
mary axes.  The  nodes  are  generated  in  rows  parallel  to  the  x axis 
and  are  numbered  consecutively  within  each  row,  in  the  direction  of  the 
increasing  x axis. 

4 = Identification  number  for  flat-plate  geometry. 

N1  = The  node  number  of  the  corner  with  the  minimum  x and 
y coordinates  and  the  maximum  z coordinate. 

NX  = The  number  of  nodes  in  the  x direction.  In  Figure  15^ 

NX  = 9. 

NY  = The  number  of  nodes  in  the  y direction.  In  Figure  15, 

NY  = 6. 

NZ  = The  number  of  nodes  in  the  z direction.  In  Figure  15, 

NZ  = 5. 

IDY  - The  node  number  increment  between  adjacent  nodes  when 
moving  parallel  to  the  y axis  in  the  positive  direction.  If 
the  entire  plate  is  formed  with  one  shape,  then  IDY  should 
be  equal  to  NX.  If,  however,  the  entire  plate  is  to  be 
formed  with  multiple  plate  shapes,  then  IDY  should  be 
equal  to  the  total  number  of  nodes  along  the  x direction  of 
the  system  of  individual  plate  shapes.  It  is  only  necessary 
to  use  multiple  shapes  when  the  expansion  factors  must 
change  along  any  of  the  three  axes. 
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Figure  15,  Geometry  of  Special  Plate  Shapes 


IDZ  = The  node  number  increment  between  adjacent  nodes 

when  moving  parallel  to  the  z axis  in  the  negative 
direction.  If  the  entire  plate  is  formed  with  one 
shape,  then  IDZ  should  be  equal  to  NX  ' NY.  For 
multiple  shapes  it  should  be  equal  to  the  total  number 
of  nodes  in  each  x-y  plane  for  the  system  of  individual 
plate  shapes. 

lY  =1  gives  restraint  in  the  y direction  if  y = 0. 

X-EXPAND  = Factor  by  which  the  distance  between  adjacent  nodes 

changes  in  the  increasing  x direction  from  XI  to  XN. 
Same  as  described  for  the  Projectile  Node  Data  Cards. 

Y-EXPAND  = Factor  by  which  the  distance  between  adjacent  nodes 

change  in  the  increasing  y direction  from  Y1  to  YN. 

Z-EXPAND  = Factor  by  which  the  distance  between  adjacent  nodes 
changes  in  the  decreasing  z direction  from  Z1  to  ZN. 

Xl  = The  X coordinate  of  node  Nl.  This  must  be  the  mini- 

mum X coordinate  of  the  plate  shape. 

Y1  " The  y coordinate  of  node  Nl.  This  must  be  the  mini- 

mum y coordinate  of  the  plate  shape. 

Z1  = The  z coordinate  of  node  Nl.  This  must  be  the  maxi- 

mum z coordinate  of  the  plate  shape. 

XN  = The  X coordinate  of  NN.  This  must  be  the  maximum 

X coordinate  of  the  plate  shape. 

YN  = The  y coordinate  of  node  NN.  This  must  be  the  maxi- 

mum y coordinate  of  the  plate  shape. 

ZN  = The  z coordinate  of  node  NN.  This  must  be  the  mini- 

mum z coordinate  of  the  plate  shape. 
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• Element  Card  for  Each  Rod  Shape  (1015)  - This  card  generates  ele- 
ments for  the  rod  shape  illustrated  in  Figure  13  and  described  by  the 
. 3 Node  Cards  for  Each  Rod  Shape. 

1 = Identification  number  for  rod  geometry, 

N1  = The  node  number  of  the  centerline  node  on  the  top  end  of 

the  rod. 

NRING  = The  number  of  rings  (1,  2,  3,  4 or  5)  in  the  rod. 

NLAY  = The  number  of  layers  of  elements  in  the  rod, 

ELEl  = The  element  number  of  the  first  element  in  the  rod.  The 

elements  are  numbered  consecutively  and  are  generated 
in  columns  of  composite  brick  elements  beginning  with 
top  layer  1 and  ending  at  bottom  layer  NLAY.  The  entire 
first  ring  of  elements  is  generated  before  the  second 
ring,  etc.  , and  the  composite  brick  elements  are  gener- 
ated in  a counterclockwise  manner  for  each  ring  of  ele- 
ments. The  total  number  of  elements  in  a rod  shape  is 
dependent  on  the  number  of  layers  (NLAY)  and  the  number 
of  elements  per  layer  (EPL)  shown  in  Figure  13. 

(Example:  If  NLAY  = 10  and  NRING  = 3,  the  total  number 
of  elements  is  NLAY  * EPL  = 10  • 96  = 960.  ) 

Ml  - M5  = Material  numbers  (1,  2^  3,  4 or  5)  for  the  rings  of  the 
rod. 


• Element  Card  for  Each  Nose  Shape  (315,  5X,  615)  - This  card  generates 
elements  for  the  nose  shapes  shown  in  Figure  14  and  described  by  the 
3 Node  Cards  for  Each  Nose  Shape. 

INOSE  = 2 identifies  a rounded  nose  shape. 

INOSE  = 3 identifies  a conical  nose  shape, 


67 


N1  = The  node  number  of  the  centerline  node  at  the  top  of  the 

nose  at  the  rod  interface. 

NRING  = Number  of  rings  in  the  nose  and  the  rod. 

ELEl  = The  element  number  of  the  first  element  in  the  nose  shape. 

The  elements  are  numbered  consecutively  from  the  inner 
ring  to  the  outer  ring. 

Ml  - M5  = Material  numbers  (1,  2,  3,  4 or  5)  for  the  rings  of  the 
nose  shape. 

• Element  Card  for  Each  Flat-Plate  Shape  (315,  5X,  515)  - This  card 

generates  elements  for  the  flat-plate  shape  illustrated  in  Figure  15  and 

described  by  the  3 Node  Cards  for  Each  Flat  Plate  Shape. 

4 = Identification  number  for  flat-plate  geometry. 

N1  = The  node  number  of  the  corner  node  shown  in  Figure  15. 

4 = Identification  number  for  flat-plate  geometry. 

ELEl  = The  element  number  of  the  first  element  in  the  flat  plate. 
The  elements  are  numbered  consecutively  and  are  gener- 
ated in  columns  of  composite  brick  elements  beginning 
at  node  Nl,  The  columns  of  elements  go  in  the  direction 
of  the  increasing  x axis.  Each  column  of  elements  is  to 
the  positive  y direction  of  the  preceding  column.  After  a 
layer  is  complete,  the  next  lower  layer  is  generated  in  a 
similar  manner. 

MATE  = Material  number  (1,  2,  3,  4 or  5)  of  the  flat  plate, 

NLX  = Number  of  layers  of  composite  brick  elements  in  the 

X direction.  The  total  number  of  nodes  along  the  x direc- 
tion must  be  NLX  + 1.  In  Figure  15,  NLX  = 8. 

NLY  = Number  of  layers  of  composite  brick  elements  in  the 

y direction.  The  total  number  of  nodes  along  the  y direc- 
tion must  be  NLY  +1.  In  Figure  15,  NLY  = 5. 
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NLZ 


Number  of  layers  of  composite  brick  elements  in  the 
z direction.  In  Figure  15,  NLZ  = 4, 


4.  2 INPUT  DATA  FOR  A RESTART  RUN 


A restart  run  is  one  which  reads  data  from  Tape2,  These  data  must  have 
been  previously  generated  by  setting  ISAVE  = 1 on  a Data  Output  Card  for 
an  initial  geometry  run  or  a previous  restart  run.  The  descriptions  which 
follow  are  for  the  input  data  in  Figure  10. 


• Identification  Card  (215.  FlO,  0)  - 


CASE  = Case  number  for  run  identification.  This  does  not  have  to 
be  identical  to  that  of  the  previous  run. 


CYCLE  = The  cycle  number  at  which  the  restart  occurs.  The  cycle 
numbers  for  which  restart  files  are  written,  are  given  in 
the  printed  output  of  the  previous  run, 

CPMAX  = Central  processor  time  (hours)  at  which  the  results  will  be 
written  onto  a restart  tape  and  the  run  will  stop.  This 
feature  will  be  bypassed  if  CPMAX  = 0. 


• Integration  Time  Increment  Card  (3F10.  0)  - 


DTMIN 


STF 


TMAX 


The  minimum  integration  time  increment  which  is  allowed. 
If  At  from  Equation  (42)  is  less  than  DTMIN,  the  results 
will  be  written  into  the  restart  tape  and  the  run  will  stop. 

The  fraction  of  the  sound  speed  transit  time  used  for  the 
integration  time  increment.  This  is  identical  to  C^  in 
Equation  (42)  and  must  be  less  than  unity. 

The  maximum  time  the  problem  is  allowed  to  run.  This 
time  refers  to  the  dynamic  response  of  the  system,  not 
the  central  processor  time  (CPMAX)  defined  in  the 
Identification  Card. 
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Data  Output  Cards  (2F10.  0,415)  - Identical  to  the  Data  Output  Cards 
for  the  Initial  Run  Input  Data  described  in  4.  1. 

4,  3 INPUT  DATA  FOR  PLOTTING  PROGRAM 

A Calcomp  plotting  program  has  been  written  to  plot  the  cross-sectional 
geometry  of  the  x-z  plane  at  y=0.  Plots  can  be  obtained  for  various  times 
by  setting  IPLOT  = 1 on  the  Data  Output  Cards.  This  will  write  the  required 
plotting  data  on  Tape3  so  they  can  be  read  by  the  plotting  program.  The 
vertical  z axis  is  currently  10  inches  long,  and  the  horizontal  x axis  is  as 
specified.  Each  plot  requires  one  card,  and  the  cards  must  be  arranged 
in  order  of  increasing  cycle  number.  The  last  card,  which  will  stop  the 
program,  should  have  four  nines  (9999)  in  columns  2-5. 

• Plot  Card  (215.  4F1Q.  0.  3A6.  IX.  II.  4X.  14. 4X.  14)  - 

CASE  = The  same  case  number  as  used  on  the  Identification  Card 
for  the  run  being  plotted.  The  last  card  in  the  data  deck 
should  have  CASE  = 9999.  (Columns  1-5) 

CYCLE  = The  cycle  number  of  the  plot  which  is  desired.  The  cycle 
numbers  of  the  data  written  on  Tape3  are  given  in  the 
printed  output  of  the  main  program,  (Column  6-10) 

ZMAX  = The  maximum  z coordinate  to  be  included  in  the  plot. 

Any  elements  having  a z coordinate  greater  than  ZMAX 
will  not  be  plotted,  (Columns  11-20) 

ZMIN  = The  minimum  z coordinate  to  be  included  in  the  plot.  Any 
elements  having  a z coordinate  less  than  ZMIN  will  not  be 
plotted.  The  vertical  z axis  is  10  inches  long  and  extends 
from  ZMAX  down  to  ZMIN.  (Columns  21-30) 

XMAX  = The  maximum  x coordinate  to  be  included  in  the  plot.  Any 
elements  having  an  x coordinate  greater  than  XMAX  will 
not  be  plotted.  (Columns  31-40) 
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XMIN 


TITLE 

IE 


IP 


IF 


The  minimiiin  x coordinate  to  be  included  in  the  plot. 

Any  elements  having  an  x coordinate  less  than  XMIN  will 
not  be  plotted.  The  length  of  the  horizontal  x axis  can 
vary.  This  axis  will  have  the  same  scale  factor  as  the 
vertical  z axis.  (Columns  41-50) 

The  title  which  is  written  on  the  plot.  The  time  and  cycle 
number  are  also  written  on  the  plot.  (Columns  51-68) 

1 will  write  "E"  at  the  center  of  all  elements  which  are 
in  the  elastic  range.  Will  not  write  if  it  is  left  blank. 
(Column  70) 

1 will  write  a at  the  center  of  all  elements  which  are 
in  the  plastic  range.  Will  not  write  if  it  is  left  blank. 
(Column  75) 

1 will  write  an  "P"  at  the  center  of  all  elements  which 
are  failed  in  shear  and  tension.  Will  not  write  if  it  is 
left  blank.  (Column  80) 


4.4  SAMPLE  INPUT  DATA 

A listing  of  input  data  for  a sample  run  is  shown  in  Figure  16.  The  projec- 
tile is  a long  rod,  with  a hemispherical  nose,  impacting  a flat-plate  target 
at  an  oblique  angle.  Consistent  units  (pound-inch-second)  are  used  for  the 
input  data.  The  projectile  (MATLl)  and  the  target  (MATL3)  are  similar 
steels,  with  the  projectile  having  slightly  higher  yield  and  ultimate  stresses. 
The  projectile  geometry  is  generated  with  one  special  rod  shape  and  one 
special  nose  shape.  The  total  length  is  4.  0 inches  and  the  diameter  is  0.  4 
inch.  The  target  requires  six  special  flat-plate  shapes  for  the  nodes  since 
the  expansion  factors  (X-EXPAND,  Y-EXPAND,  Z-EXPAND)  are  not  con- 
stant for  the  X and  y directions.  Only  one  special  element  card  is  required 
for  the  plate,  however.  The  length  of  the  plate  is  8.  0 inches,  the 
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Figure  16,  Input  Data  for  the  Sample  Problem 
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half -width  is  3.  0 inches  and  the  thickness  is  0.4  inch.  This  problem  con- 
tains 1992  nodes  and  8304  elements.  The  master  sliding  surface  is  defined 
by  five  rows  of  21  nodes  for  a total  of  105  master  nodes.  A total  of  108 
surface  nodes  on  the  nose  and  lower  portion  of  the  rod  are  designated  as 
slave  nodes.  The  impact  conditions  are  60  degrees  from  normal  at  a net 
velocity  of  60,  000  inches/ second,  A Calcomp  plot  is  shown  in  Figure  17 
for  the  cross-sectional  geometry  at  cycle  1, 


4.  5 OUTPUT  DATA  DESCRIPTION 

The  output  data  are  generally  clearly  described  using  the  terminology  of 
this  report,  A summary  of  the  printed  output  for  an  initial  run  is  as  follows 

• Material  data  are  printed  for  the  input  and  generated  data. 

• All  the  geometric  input  data  are  printed  in  a form  similar 
to  that  used  to  read  the  data. 

• For  each  element,  the  four  nodes,  the  volume  and  material 
number  are  printed. 

• For  each  node,  the  initial  coordinates,  mass  and  restraints 
are  printed.  The  slave  nodes  are  identified  by  adding  a "l" 
before  the  XYZ  restraint.  (Example:  If  the  XYZ  restraint 
is  110,  the  node  is  restrained  in  the  x and  y directions.  If 
the  XYZ  restraint  is  1110,  it  is  a slave  node  with  the  same 
restraints.  ) 

• A summary  of  the  initial  geometry  and  impact  conditions  is 
printed.  Included  are  the  number  of  nodes  and  elements, 
mass,  eg  positions,  velocities,  momenta,  kinetic  energies 
and  integration  time  increment  data. 
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• For  each  integration  time  increment  selected  data  are  printed. 
These  consist  of  cycle  number  (CYCLE),  the  time  (TIME),  the 
integration  time  increment  (DT),  the  element  number  which 
governs  the  time  increment  (ECRIT),  the  total  system  kinetic 
energy  (TOTAL  K.  E.  ),  the  projectile  momenta  (PROJ  MVX, 

PROJ  MVY,  PROJ  MVZ),  the  plastic  work  in  the  projectile 
(PROJ  PLAST),  the  plastic  work  in  the  total  system  (TOT  PLAST), 
and  the  total  energy  in  the  system  (TOT  ENERGY).  If  ECRIT  = 0, 
the  integration  time  increment  calculated  by  Equation  (42)  is 
either  greater  than  DTIMAX,  or  greater  than  1.  1 times  the  inte- 
gration time  increment  for  the  previous  cycle.  If  the  total 
energy  in  the  system  remains  constant,  the  approximate  formu- 
lation for  the  internal  energy  in  Equation  (56)  is  adequate. 

• A summary  of  data  for  the  projectile,  target  and  total  system 
(projectile  and  target)  is  printed  for  selected  times  specified  on 
the  Data  Output  Cards.  These  data  consist  of  mass,  eg  positions, 
linear  momenta,  linear  velocities,  angular  momenta  and  angular 
velocities. 

• Node  data  are  printed  at  the  selected  times  specified  on  the  Data 
Output  Cards  if  ECHECK  is  less  than  1000.  These  data  include 
the  node  number  (NODE),  the  coordinates  (X,  Y,  Z),  the  velocities 
(XDOT,  YDOT,  ZDOT)  and  the  accelerations  (XDD,  YDD,  ZDD). 

If  the  acceleration  data  are  not  printed,  this  indicates  that  the 
node  is  a slave  node  which  is  currently  in  contact  with  the  master 
surface. 

• Element  data  may  be  printed  at  the  selected  times  specified  on 
the  Data  Output  Cards  if  ECHECK  is  less  than  1000.  If  ECHECK 
is  negative,  only  the  plastic  and  failed  element  data  are  printed. 
For  positive  ECHECK  less  than  1000,  only  those  elements  with 
an  equivalent  strain  equal  to  or  greater  than  ECHECK,  will  have 


their  data  printed.  If  the  element  is  plastic  range,  these  data 
consist  of  the  element  number  (ELE),  the  normal  strains 
(EX,  EY,  EZ)  and  shear  strains  (EXY,  EXZ,  EYZ)  given  by 
Equations  (7)  to  (12),  the  equivalent  strain  (EBAR)  given  by 
Equation  (14),  the  equivalent  strain  rate  (EDOT),  the  volume- 
tric strain  (DVOL)  given  by  Equation  (13),  the  equivalent  tensile 
strength  (SEFF)  given  by  Equation  (31),  the  hydrostatic  pressure 
(PRESURE)  given  by  Equation  (32),  the  artifical  viscosity 
[Q(VIS)]  given  by  Equation  (33),  the  normal  de viator  stresses 
(SX,  SY,  SZ)  given  by  Equations  (25)  to  (27),  and  the  shear 
stresses  (SXY,  SXZ,  SYZ)  given  by  Equations  (28)  to  (30), 

If  the  element  is  in  the  elastic  range,  EDOT  is  set  equal  to  1, 
and  SEFF  is  the  equivalent  stress  given  by  Equation  (21). 

Also,  the  normal  elastic  stresses  of  Equations  (15)  to  (17)  are 
presented  in  the  form  of  a hydrostatic  pressure  and  deviator 
stresses. 

• Internal  loads  can  be  obtained  for  slender  projectiles  at  selected 
times  if  the  nodal  geometry  is  consistent  with  the  geometry 
generator  for  rods  and  if  ILOAD  = 1 on  the  Data  Output  Cards. 
The  output  is  given  along  the  centerline  of  the  projectile  from 
nodes  NTOP  to  NBOT  as  specified  on  the  Slide/  Load  Card. 
Included  are  the  node  numbers  (NODE),  the  coordinates  of  the 
nodes  (X,  Z),  the  bending  moment  (BENDING  MOMENT),  the  x 
and  z forces  (X  LOAD,  Z LOAD)  the  transverse  shear  force 
(SHEAR  LOAD)  and  the  axial  force  (AXIAL  LOAD).  The  bending 
moment  is  calculated  at  the  node  positions  and  acts  on  the  lower 
(forward)  end  of  the  upper  (aft)  portion  of  the  projectile  as  shown 
in  Figure  18.  A clockwise  moment  is  positive  when  looking  along 
the  y axis.  The  x and  y forces  are  calculated  for  the  increments 
between  the  nodes  and  are  positive  when  applied  in  the  positive 
x and  z directions.  The  shear  force  is  positive  when  it  acts  in  a 
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Figure  18.  Sign  Convention  for  Internal  Projectile  Loads 


clockwise  direction  about  the  upper  (aft)  end  of  the  projectile 
when  looking  along  the  y axis.  The  axial  force  is  positive  when 
the  projectile  is  in  tension.  The  results  are  valid  only  along 
the  free  portion  of  the  projectile,  and  are  not  valid  for  the  por- 
tions which  are  in  contact  with  the  target. 

4.  6 DIAGNOSTICS 

The  following  diagonistics  are  provided. 

• CP  TIME  LIMIT  EXCEEDED  - This  means  the  central  processor 
time  is  greater  than  specified  by  CPMAX,  on  the  Identification 
Card.  The  results  are  written  on  the  restart  tape  and  the  run 

is  discontinued. 

• MINIMUM  TIME  INCREMENT  HAS  BEEN  VIOLATED  - This 
means  the  time  increment  calculated  by  Equation  (42)  is  less 
than  specified  by  DTMIN  on  the  Integration  Time  Increment 
Card,  The  results  are  written  on  the  restart  tape  and  the 
run  is  discontinued.  If  it  is  desired  to  continue  the  run  to 
later  times,  DTMIN  must  be  decreased  for  the  next  restart 
run. 
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• ENERGY  CHECK  INDICATES  NUMERICAL  INSTABILITY  - 
This  means  the  kinetic  energy  of  the  system  is  at  least  10  per- 
cent greater  than  the  initial  kinetic  energy.  Since  this  is  not 
physically  possible,  the  run  must  be  numerically  unstable. 

• END  OF  RESTART  TAPE  - This  means  a restart  run  has 
been  attempted,  but  the  restart  tape  does  not  have  data  for 
the  cycle  number  specified  on  the  Identification  Card, 

4,  7 CENTRAL  PROCESSOR  TIME  ESTIMATES 

For  a limited  number  of  problems  run  on  a Honeywell  6080  Computer,  the 
program  has  used  approximately  0,  025  central  processor  second  per  node 
per  cycle. 


4.  8 CENTRAL  MEMORY  STORAGE  REQUIREMENTS 
AND  ALTERATIONS 

The  program  is  currently  limited  to  five  materials,  2000  nodes  and  two 
sliding  surfaces,  each  with  a maximum  of  300  slave  nodes  and  300  master 
nodes.  The  element  blocks  contain  data  for  200  elements.  This  requires 
69K  Central  Memory  Storage  on  a Honeywell  6080  Computer,  To  increase 
the  capability  of  the  program  for  problems  involving  more  than  2000  nodes, 
it  is  necessary  to  increase  the  size  of  the  node  arrays  in  Subroutines  ELEG, 
GEOM,  GPLOT,  LOAD,  LOOP,  NODEG,  RECALL,  SAVE,  SDATA,  SLIDE 
and  START,  The  sliding  sxarface  arrays  are  contained  in  Subroutines 
GEOM,  LOOP,  RECALL,  SAVE  and  SLIDE,  If  it  is  desired  to  change  the 
block  size  for  the  elements,  the  element  arrays  are  contained  in  Subroutines 
ELEG,  GEOM,  GPLOT,  LOOP,  RECALL,  SAVE  and  SLIDE.  If  the  size 
of  the  element  arrays  is  changed,  it  is  necessary  to  redefine  NBSIZE  (Card 
37  in  Subroutine  GEOM)  to  be  equal  to  the  dimension  of  the  element  arrays. 
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